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Abstract 

We introduce a new particle-based hybrid code for planetary accretion. The 
code uses an A-body routine for interactions with planetary embryos while 
it can handle a large number of planetesimals using a super-particle approx¬ 
imation, in which a large number of small planetesimals are represented by a 
small number of tracers. Tracer-tracer interactions are handled by a statis¬ 
tical routine which uses the phase-averaged stirring and collision rates. We 
compare hybrid simulations with analytic predictions and pure A-body sim¬ 
ulations for various problems in detail and hnd good agreements for all cases. 
The computational load on the portion of the statistical routine is compara¬ 
ble to or less than that for the A-body routine. The present code includes 
an option of hit-and-run bouncing but not fragmentation, which remains for 
future work. 

Keywords: Accretion; Planetary formation; Planetary Rings; Planets, 
migration; Origin, Solar System 


1. Introduction 

Terrestrial planets and cores of giant planets are generally considered to 
have formed through accretion of many small bodies called planetesimals. 
To simulate accretion processes of planets, two methods, which are comple¬ 
mentary to each other, have been applied. 

The hrst one is A-body simulations in which orbits of all bodies are 
numerically integrated and gravitational accelerations due to other bodies 
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are calculated in every time step (e.g., Kokubo and Ida, 1996; Chambers 
and Wetherill, 1998; Richardson et al., 2000; Morishima et ah, 2010). N- 
body simulations are accurate and can automatically handle any complicated 
phenomena, such as resonant interactions and spatially non-uniform distri¬ 
butions of planetesimals. Gravity calculations are accelerated by such as 
tree-methods (Richardson et ah, 2000) or special hardwares (Kokubo and 
Ida, 1996; Grimm and Stadel, 2014), and a large time step can be used 
with sophisticated integrators, such as Mixed Variable Symplectic (MVS) or 
Wisdom-Holman integrators (Duncan et ah, 1998; Chambers, 1999). Even 
with these novel techniques, iV-body simulations are computationally intense 
and the number of particles or the number of time steps in a simulation is 
severely limited. 

The second approach is statistical calculations in which planetesimals are 
placed in two dimensional (distance and mass) Eulerian grids, and the time 
evolutions of the number and the mean velocity of an ensemble of plan¬ 
etesimals in each grid are calculated using the phase-averaged collision and 
stirring rates (e.g., Greenberg et al., 1978; Wetherill and Stewart, 1989, 1993; 
Inaba et ah, 2001; Morbidelli et ah, 2009; Kobayashi et ah, 2010). While 
this approach does not directly follow orbits of individual particles, it can 
handle many particles, even numerous collisional fragments. Largest bodies, 
called planetary embryos, are handled in a different manner than small bod¬ 
ies, taking into account their orbital isolation. The mutual orbital separation 
between neighboring embryos is usually assumed to be 10 mutual Hill radii. 

The last assumption is not always guaranteed, particularly in the late 
stage of planetary accretion (e.g., Ghambers and Wetherill, 1998). To handle 
orbital evolution of embryos more accurately, Eulerian hybrid codes 0 have 
been developed (Spaute et ah, 1991; Weidenschilling et al., 1997; Bromley 
and Kenyon, 2006; Glaschke et al., 2014), in which small planetesimals are 
still handled by the Eulerian approach whereas orbital evolutions of largest 
embryos are individually followed by such as iV-body integrations. Gravita¬ 
tional and collisional interactions between embryos and small planetesimals 
are handled using analytic prescriptions. Glaschke et al. (2014) took into 
account radial diffusion of planetesimals due to gravitational scattering of 


^The term ’’hybrid” means a combination of two different approaches/schemes. While 
this term might be used as a combination of Eulerian grids and Lagrangian embryos by 
the authors, we use the term for a combination of statistical calculations and Wbody 
calculations. 
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embryos and their code can approximately handle gap opening around em¬ 
bryos’ orbits. 

A Lagrangian hybrid method has also been introduced by Levison et ah 
(2012) (LDT12 hereafter). In their LIPAD code, a large number of planetes- 
imals are represented by a small number of Lagrangian tracers. This type 
of approach is called a super-particle approximation and is also employed in 
modeling of debris disks (Krai et ah, 2013; Nesvold et ah, 2013) and planetes- 
imal formation (Johansen et ah, 2007; Rein et ah, 2010). Orbits of individ¬ 
ual tracers are directly followed by numerical integrations, and interactions 
between planetesimals (stirring and collisions) in tracers are handled by a 
statistical routine. Embryos are represented by single particles and the ac¬ 
celerations of any bodies due to gravity of embryos are handled in the iV-body 
routine. Lagrangian hybrid methods have several advantages than Eulerian 
hybrid methods. For example, Lagrangian methods can accurately handle 
spatial inhomogeneity in a planetesimal disk (e.g., gap opening by embryos), 
planetesimal-driven migration, resonant interactions between embryos and 
small planetesimals, and eccentric ringlets. A drawback of Lagrangian meth¬ 
ods would be computational cost as orbits of all tracers need to be directly 
integrated. Therefore, it is desirable that Lagrangian hybrid methods can 
handle various problems accurately even with a moderately small number of 
tracers. 

In this paper, we develop a new Lagrangian hybrid code for planet for¬ 
mation. While we follow the basic concept of particle classes introduced by 
LDT12, recipes used in our code are different from those used in the LIPAD 
code in many places. The largest difference appears in the methods to han¬ 
dle viscous stirring and dynamical friction. The LIPAD code solves pair-wise 
interactions while our method gives the accelerations of tracers using the 
phase-averaged stirring and dynamical friction rates. While the LIPAD code 
conducts a series of three body integrations in the stirring routine (in the 
shear dominant regime) and in the routine of sub-embryo migration during 
a simulation, our code avoids them by using unique routines. The complete 
list of comparison with the LIPAD code turns out to be rather long and is 
given in Appenedix G. 

In Section 2, we explain our method. In Section 3, we show various tests 
of the new hybrid code and compare them with analytic estimates and pure 
iV-body simulations. The computational speed and limitations of our code 
are discussed in Section 4. The summary is given in Section 5. For the sake 
of clarity, specihc derivations are deferred to appendices. 
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Fig. 1. Schematic illustration of our hybrid code. Accelerations due to gravity 
of embryos are handled by the A^-body routine. Tracer-tracer interactions 
and back reaction of tracers on sub-embryos including collision are handled 
by the statistical routine. 

2. Method 

2.1. Particle classes 

The particle classes in our hybrid code are the same as those introduced 
by LDT12. The code has three classes of particles: tracers, sub-embryos, and 
full-embryos (Fig. 1). In the present paper, we do not consider fragmentation 
of planetesimals and dust-tracers are not introduced. A tracer represents a 
large number of equal-mass planetesimals on roughly the same orbits. The 
mass of a planetesimal and the number of planetesimals in the tracer i {i 
for the index of the tracer) are dehned to be rrii and /cj. Therefore, the 
tracer mass is given by kifrii. Through collisional growth, m* increases and 
ki decreases while kifrii remains close to its original value. We allow mass 
exchanges between tracers through collisions so the tracer mass, kifrii, is 
not necessarily hxed but has the upper and lower limits, rut^max and 

< kifrii < mi,max)- We employ rrii^max = 2.0mio and = O.lmto 
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in this paper, where rrito is the minimum mass of a sub-embryo. The mass 
rrito is usually the initial mass of a tracer, which is the same for all tracers at 
the beginning of a simulation. If ki = 1 and m* > rrito, the tracer is promoted 
to a sub-embryo. If mt > ffmto, the sub-embryo is further promoted to a 
full-embryo, where we employ the numerical factor /f of 100, as recommended 
by LDT12. The number ki is an integer in our model. 

Orbits of any types of particles are directly integrated. We use a Mixed 
Variable Symplectic integrator known as SyMBA (Duncan et ah, 1998), 
which can handle close encounters between particles. This integrator is also 
used by LDT12. The collisional and gravitational interactions between full- 
embryos with tracers are directly handled in every time step of orbital in¬ 
tegrations, as is the case of pure V-body simulations. On the other hand, 
tracer-tracer interactions are handled in a statistical routine which is de¬ 
scribed in subsequent sections in great detail. While the time step for the 
orbital integration 6t is ~ 10“^ of the orbital period, the time step At for 
the statistical routine can be taken to be much larger as long as At is suf- 
hciently smaller than the stirring and collisional timescales. In this paper, 
we employ At = 306t for all test simulations. We conhrmed that almost the 
same outcome is obtained even with a smaller At for all the test simulations. 

In interactions between a tracer and a full-embryo in the V-body routine, 
the tracer is assumed to be a single particle with the mass equivalent to the 
total mass of planetesimals in the tracer. This is not a problem as long as 
the embryo is sufficiently massive compared with tracers. However, if the 
embryo mass is similar to tracer masses as is the case immediately after its 
promotion, and if embryo-tracer interactions are handled by the direct N- 
body routine, the embryo suffers artihcially strong kicks from tracers. To 
avoid this issue, sub-embryos are introduced. Accelerations of sub-embryos 
due to gravitational interactions with tracers are handled by the statistical 
routine whereas accelerations of tracers due to gravitational interactions with 
sub-embryos are handled by the direct V-body routine (Fig. 1). Collisions of 
planetesimals in tracers with sub-embryos are also handled in the statistical 
routine to avoid artihcially large mass jumps in sub-embryos, contrary to 
LDT12 (see Appendix G for more discussion). 

2.2. Neighbor search for calculating surface number densities 

Our statistical routine uses the surface number density of nearby tracers. 
We hrst describe how we determine the surface number density for tracer- 
tracer interactions. The surface number density of tracers for interactions 
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Fig. 2. Schematic illustration of neighboring tracer search. The hlled circle 
is the target tracer and open circles are interloping tracers (interlopers). The 
surface number density of planetesimals in an interloper is simply given by 
the number of planetesimal in the interloper divided by the area of the region 
between two arcs. 


between a sub-embryo and tracers is derived in a similar way, but slightly 
modified so that close and distant encounters are treated differently (Sec¬ 
tion 2.5). 

Consider a target tracer i, surrounded by nearby interloping tracers. The 
cylindrical co-ordinate system (r, 9, z) is introduced and we consider a curved 
region, called the region i, located at z = 0 and centered at the position of 
the tracer i (Fig. 2). The radial half width and the angular half width in the 
6 direction of the region i are given by Sr and S6, respectively. The area of 
this region is given by 

Si = AviSrSO. ( 1 ) 


For neighboring tracer search, tracers are sorted into 2D (r, 6) grids and only 
nearby cells of the target tracer are checked whether other tracers are in the 
designated region. If the interloper j is in the region i, the surface number 
density of planetesimals in the interloper j around the target i is simply 
expressed as 



( 2 ) 
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The density Uj is an instantaneous density. At other times, the interloper j 
may not be in the region i but other interlopers having planetesimal masses 
similar to mj may be in the region i. After averaging over time, the encounter 
rate with planetesimals in a certain mass range is expected to converge re¬ 
gardless of the size of Si. 

A small 6r is favorable to capture a small radial structure accurately. 
However, it is not reasonable to use 6r smaller than a typical radial length 
of gravitational influence of a planetesimal, which is usually about ~ 

Here, Ru is the Hill radius given by 


Rh — O^i 



( 3 ) 


where ai is the semi-major axis of the tracer i and Mq is the mass of the 
central star. Since the mass of a possible largest planetesimal is rrito, we 
employ 

for tracer-tracer interactions. For interactions between a sub-embryo and 
tracers, a larger Sr needs to be used and we employ 10 Hill radii of a possible 
largest sub-embryo with the mass of ftmto- 

An advantage of the choice of Sr given by Eq. (II]) is that we can omit 
the procedure for orbital isolation of planetesimals. Orbital isolation of the 
largest bodies occurs if the sum of orbital spaces of these bodies is less than 
the width of the region of interest (Wetherill and Stewart, 1993). If a plan¬ 
etesimal is assumed to occupy an orbital space of 10 Hill radii (Kokubo and 
Ida, 1998) and if the tracer masses are m^o, the sum of the orbital spaces 
of planetesimals in a pair of tracers are always larger than 2Sr. Since our 
model allows a tracer mass lower than rrito, isolation can potentially occur if 
the tracer masses and the numbers of planetesimals in the tracers are both 
small. Such cases are rather rare and we can ignore them. Therefore, or¬ 
bital isolation can occur only for embryos and that is handled by the A^-body 
routine. An additional discussion is given in Appendix G. 

A small SO is favorable to handle non-axisymmetric structures, in which 
the surface density is not azimuthally uniform. Such structures usually ap¬ 
pear if external massive bodies exist (Kortenkamp and Wetherill, 2000; Na- 
gasawa et ah, 2005; Levison and Morbidelli, 2007; Queck et ah, 2007). A 
small S6 can also be chosen, if we want to suppress the number of interlopers 
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to save computational costs (Section 4.1). In this paper, we usually adopt 
69 = O.Svr, but dependency on 69 is also examined. 

Even if the tracer j is in the region i, the tracer i is occasionally outside 
the region j used in neighboring search around the tracer j. To make sure 
that the tracer i and the tracer j are mutually counted (or excluded) as 
interlopers, the number density nj is evaluated using Eq. ([2]) only if ruj > rrij 
or i < j for ruj = rrij. Otherwise and only if the tracer i is in the region j, 
rij is given by kj / Sj , where Sj is the area of the region j. 

Although we only considered encounters with interlopers, planetesimals 
in the target i may encounter with other planetesimals in the same target if 
ki> 1. Appendix A describes how we correct this effect. This effect is found 
to be too small to be identihed. 

2.3. Hill’s approximation and quasi ergodic hypothesis 

Let the semimajor axis, the orbital eccentricity, the orbital inclination, 
the longitudes of pericenter and ascending node, and the mean longitude of 
the target be Oj, Cj, ij, cuj, fij, and Aj, and those for the interloper be Oj, 
ej, ij, zuj, flj, and Aj. The difference between the semimajor axes is dehned 
as dij = Oj — Oi. When the target encounters with the interloper, they may 
collide with each other at a certain probability. Their orbital elements are 
also modihed by mutual gravitational scattering. To handle these processes 
in the statistical routine in a simple manner, we make following assumptions. 

1. We use the collision probability and the mean change rates of orbital 
elements derived by previous studies which usually solved Hill’s equa¬ 
tions. The equations of motion are reduced to Hill’s equations if Hill’s 
approximation is applied (Hill, 1878; Nakazawa et ah, 1989). The con¬ 
ditions to apply Hill’s approximation are 

(a) Cj, ij, Cj, ij 1, 

(b) ruj, ruj <C M©, and 

(c) \dij\^ai. 

2. The semimajor axis differences dij and the phases of the interlopers, 
cuj, Hj, Aj are uniformly distributed. Namely, if a certain interloper in 
the region i is found, we assume that other interlopers with the same 
Cj, Zj, but different dij, zUj, flj, and Aj are distributed uniformly over 
the region i. 

3. The longitudes, zui and Hj, of the target change randomly and large 
enough on the evolution timescales of e* and Zj. This approximation is 


unnecessary if Approximation 2 is valid and either Cj or ii ij, 
because the collision and stirring rates do not depend on zuj and hi* in 
such a situation. 

Assumption 1 (Hill’s approximation) seems to work well in most of cases 
of planetary accretion. Greenzweig and Lissauer (1990) showed that the col¬ 
lision probability of a planet with planetesimals on circular and co-planar 
orbits agrees with that based on Hill’s approximation within 1% as long 
as rrii/MQ < Our test simulations for late stage planetary accretion 

(Section 3.5) indicate that Hill’s approximation also works well even at mod¬ 
erately high orbital eccentricities {cj ~ 0.3). 

Assumption 2 is acceptable if many planetesimals mutually interact via 
gravitational scattering and collisions without orbital isolation. These inter¬ 
actions lead to randomization of spatial positions and phases. We extensively 
adopt these uniformities even if the target is a sub-embryo isolated from other 
embryos, but with some approximated cares (Section 2.6). 

Assumption 3 is likely to be reasonable since the timescale of phase change 
is comparable to the stirring timescale in general (e.g, Tanaka and Ida, 1996). 
With Assumptions 2 and 3, we can assume that the time averaged stirring 
and collision rates of an individual target are equivalent to the stirring and 
collision rates averaged over phase space. While this equivalence resembles 
the so-called ergodic hypothesis, the difference is that we assume the equiv¬ 
alence on a short timescale. Although Assumption 3 may not be strictly 
correct, it greatly simplihes the problem and it is worthwhile to examine its 
validity or usefulness. 

Integrations over distributions of orbital elements of interlopers can be 
transformed to integrations over relative orbital elements of interlopers with 
respect to the target (Nakazawa et ah, 1989). The eccentricity and inclination 
vectors of the target i are given as 

Si = (d cos zui, Ci sin Wi), 
ii = (fiCosr2j,fiSinr2j). 

The same characters but with the subscript j represent the eccentricity and 
inclination vectors for the interloper j. The relative eccentricity and inclina¬ 
tion vectors are given as (Fig. 3) 


Gr — (^Cj-COS ZUJ-^ Cj-sin ZUy') — Gj — Gj, 
ir = {ir COS rir, V sin Gr) = ij ~ ii, 


(6) 
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Fig. 3. Definitions of the eccentricity vectors. The inclination vectors are 
dehned in a similar manner. 


where and Dr are the longitudes for the relative vectors. 

Nakazawa et al. (1989) adopted Assumptions 1-3 and derived the number 
of planetesimals AA^- colliding with the target i during a time interval At as 
[their Eq. (49)] 


(.?)) = 


TljOjj^jhj^jUJY^ I ir, 62, ) )-Pcol(6r 7 5 


where (e|) and (t|) are the mean squares of ej and ij of interlopers, rij is 
the surface number density of interlopers, aij = {ai -|- aj)/2 is the mean 
semi-major axis, ojk is the orbital frequency at ajj, /(cr, ir, CiAi) (e|), (f|)) 
is the distribution function of Cr and normalized as J fcterdir = 1, and 
-Pcoi(er,h) is the non-dimensional collision rate averaged over dij, Wr, and Dr 
(Appendix B.l). In the above, the reduced mutual Hill radius, hij, is given 


as 


hij 


f mi + rrij \ 

V 3Mo ) 


( 8 ) 
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The expression of / for the case of the Rayleigh distributions of Cj and ij 
is given by Eq. (35) of Nakazawa et al. (1989). In our method, on the other 
hand, / for a certain interloper is given by the Dirac delta function. Thus, 
the expected number of planetesimals AiVj in the interloper j colliding with 
the target i during At is given as 

(^r; h) '^j®pA-^KR^ol(CrAr)) (9) 

where Uj is given by Eq. ([2]). The stirring rates can be derived in a sim¬ 
ilar manner using the non-dimensional stirring rates (Section 2.5). In our 
method, the averaging over distributions of eccentricities and inclinations of 
interlopers will be made by averaging over multiple interlopers. Although it 
is not necessary to assume the distribution functions of eccentricities and in¬ 
clinations of interlopers in calculations of the collision and stirring rates, our 
stirring routine naturally reproduces Rayleigh distributions (Section 3.2). 

If the system is perturbed by a massive external body and the forced 
eccentricity is much larger than free eccentricities, the phases, vjj and Dj, 
may not be uniformly distributed. Even in such a case, our approach may 
still be applied by replacing the eccentricity and inclination vectors by the 
free eccentricity and inclination vectors whose directions are uniformly dis¬ 
tributed, provided that the forced eccentricity vectors are roughly the same 
between the target and interlopers. Thus, while uniformities of and Dj. are 
assumed, we do not necessarily assume uniformities of the relative longitudes 
Wij = Wj — Wi and Qij = Qj — Dj. The averaging over cijjj and Qij is done 
by averaging over interlopers [uj = nj(er, h, ^q)]- While this relaxing of 
the assumption does not alter Eq. (E]), it is important for dynamical friction 
(Section 2.5). The forced eccentricity vectors are not necessarily the same 
between the target and interlopers, for example, if the mass-dependent dis¬ 
sipative force works. In such a situation, our approach is not valid, although 
it is still worthwhile to examine how inaccurate it is. 

2-4^. Collision 

2.4-1. Judgement of collision 

The expected change rate of the mass rui of a planetesimal in the tracer 
i due to collisions with planetesimals in the interloper j is simply given by 
multiplying the planetesimal mass nij to Eq. ([9]) (Ida and Nakazawa, 1989; 
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Greenzweig and Lissauer, 1990; Inaba et al., 2001) 


= njmjO^jh^ijUKPcoh ( 10 ) 

j 

To avoid double counting, we only consider a case for > rrij or i < j for 
rrii = rrij. Using the mass change rate, the probability that the mass of a plan- 
etesimal in the target i increases by Am^ due to a collision with a planetesimal 
in the interloper j during the time step At is given as {dmi/dt)j{At/Ami). If 
this probability is larger than a random number which takes between 0 and 
1, we assume that all planetesimals in the target i collide with planetesimals 
in the interloper. The same procedure is applied to all other interlopers in 
the region i to check whether collisions occur between planetesimals in the 
target and those in interlopers. 

If a collision occurs and mutually colliding bodies merge, the changes 
in the masses and the numbers of planetesimals in the tracers are given as 
follows: 



mjft = m, + Am, and kib = ki (for target), (11) 

/\777 . 

rrijb = rrij and kjb = kj - -ki (for interloper), (12) 

rrij 

where the characters with the additional subscript b mean those after the 
merging. Note that the total mass is conserved in the collisions: kibrnu, + 
kjbrribj = ki-rrii + kjrrij. 

Before determining whether a collision occurs or not, we need to choose 
an appropriate Amj. The mass Arrii is usually the mass of the interloping 
planetesimal rrij, but can take different values for convenience as follows: 


Arrii 


Max[rrij,mt{0.01rrii/rrij)rrij] Case (a); kjb > 0, 
kjrrij/ki Case (b); kjb = 0 , 


(13) 


where int() means the decimals are dropped off for the number in the paren¬ 
theses. The form in the upper row of Eq. flTdll means that if rrij is smaller 
than the threshold mass, we do not consider individual collisions but the ac¬ 
cumulated mass change of the target planetesimal. We adopt the threshold 
mass of O.Olm,. We hrst estimate the post-collision mass of the interloper 
kjbrrij using Eq. flT^ and Arrii for Case (a) of Eq. flTSl) for any interlopers. 
If the resulting mass, kjbrrij, is found to be smaller than the lower limit of 
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the tracer mass, min, we use Am, of Case (b) given by the lower row in 
Eq. (IT^ instead of that of Case (a). If a collision occurs for Case (b), all the 
planetesimals in the interloper j merge with those in the target i and the 
interloper is deleted. Therefore, the mass of any tracer is kept to be lager 
than While we adopt = O.lmto, this value can be arbitrarily 

chosen: a small lower limit may be able to give accurate mass evolution while 
increase of the number of tracers can be suppressed with a large lower limit. 

2.4-2. Mass and orbit changes due to merging 

The procedure to change masses and orbital elements of planetesimals due 
to collisional merging is schematically summarized in Fig. 4. For Case (a) in 
Eq. (lT3l) . in which kji, > 0 and kj^mj > mt^rain, only some of planetesimals 
in the interloper j collide with planetesimals in the target i. In this case, 
we split the interloper j into two tracers: the interloper j and the impactor. 
The number of planetesimals in the interloper after the split is given by 
kjb in Eq. flT^ and the position and velocity of the post-split interloper are 
unchanged during the collision between the target and the impactor. For 
Case (b), all the planetesimal in the interloper are involved in the collision 
and we call the interloper the impactor. In both cases, the planetesimal mass 
and the number of planetesimals in the impactor are given by Aruj and fcj, 
respectively, and the collision process is almost the same except some small 
differences. 

When the target and the impactor are judged to collide, they do not 
physically contact each other. Even their stream lines do not usually cross 
each other. Since we assume uniformity in semimajor axis of an interloper, 
it is possible to make two stream lines cross each other by changing the 
semi-major axis of the impactor without changing other orbital elements. 
The crossing point of the stream lines is the location where the impact takes 
place. This impact point can be found as follows. First, 6 and z/r of the 
target and the impactor are aligned by changing the true anomaly of the 
target /* by A/j and that for the impactor fj by Afj. The angular changes, 
Afi and Afj, are derived by solving the following set of equations: 


(14) 

(15) 


+ /i + Afi — Wj -\- fj + Afj — 6*col, 
sin ii sin {9co\ — fli) = sin ij sin (^coi — fij), 


where ^coi is the longitude of the impacting position measured from the refer¬ 
ence direction 9 = 0. Eqs. ffTT|) and ffTSjl give two solutions for A/* mutually 
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New Tracer New Tracer Interloper 


Fig. 4. Schematic illustration of the procedure to handle collisional merging 
between planetesimals in two tracers (the target and the interloper). The 
figure shows how the masses of planetesimal (m^ and rrij), the numbers of 
planetesimals in the tracers (/c* and kj), the semimajor axes (oj and aj), and 
the orbital eccentricities (e^ and ej) change during the procedure. 
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shifted by vr, and we choose one with a smaller absolute value. Then, the true 
anomalies are changed by moving the colliding pair along their Keplerian or¬ 
bits. Second, the semimajor axis of the impactor is scaled by Vi/rj. This is 
done by scaling the impactor position by Vi/rj and the velocity by \/rJJri- 
The positions of the target and the impactor become almost identical through 
these procedures without changing their eccentricity and inclination vectors. 

Then, the merging between the target and the impactor takes place. The 
position and velocity of the target after the merging with the impactor are 
given by the mass weighted mean values (i.e., the mass center for the po¬ 
sition). After that, the semi-major axis of the target is adjusted to a^b by 
scaling the position and the velocity so that the z-component of the total or¬ 
bital angular momentum of the target and the interloper is conserved during 
the merging: 


- e|) cos in, = 

kimi\J aj(l — ef) cosi, + {kj — kj},)mjsjaj{l — e^) cosij, 

where en and in are the orbital eccentricity and inclination of the target after 
the merging. 

If the target tracer mass is equal to or exceeds the upper limit, mt,max(= 
2mio), we split it into two tracers. If kn = 2 and mn > 2mi^min, we also 
split the tracer to avoid simultaneous emerging of two sub-embryos. The 
numbers of planetesimals in the split tracers are given by kn/‘2, if kn is even, 
and kn/‘2 ± 0.5 if kn is odd. We give weak velocity changes, symmetrically 
in momenta, to the two new tracers so that their relative eccentricity and 
inclination are about hij. This will avoid an encounter velocity artihcially 
too low between them in the next encounter. The semi-major axes of the 
two tracers are adjusted so that the orbital separation between two tracers 
is randomly chosen between 0 to 6r but the ^-component of the total orbital 
angular momentum is conserved. For Case (b), we set the semi-major axis 
of one of the new tracers be the one for the deleted interloper instead of 
choosing the orbital separation randomly. Finally, the true anomalies of the 
two tracers (or one if no splitting is made) are changed by — A/j and — A/j. 

Through the processes described above, tracer masses remain between 
rrit^rm-a and mt^niax- The number of tracers remains unchanged or increases for 
Case (a) of Eq. flT^ whereas it remains unchanged or decreases for Case (b). 
As a result, the total number of tracers remains similar to the initial value 
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even after many collisions between planetesimals. Our several experiments 
show that increase of the tracer number from the initial value is suppressed 
less than 30% in most cases. Our current code also has an option of hit-and- 
run collisions and the procedure is described in detail in Appendix E. In this 
case, neither the number of planetesimals nor that of tracers changes. 


2.5. Stirring and radial diffusion 

2.5.1. Rates of stirring and radial diffusion 

The orbital elements of planetesimals evolve due to gravitational inter¬ 
actions between them. The interactions are separated into viscous stirring 
(VS) and dynamical friction (DF). Viscous stirring increases eccentricities 
and inclinations in expense of the tidal potential (i.e., radial diffusion) and 
dynamical friction leads to energy equipartition between particles (Ida, 1990, 
Stewart and Ida, 2000; Ohtsuki et ah, 2002). Dynamical friction can be sep¬ 
arated into stirring and damping parts (DFS and DFD) and sometimes the 
latter part only is called dynamical friction (Binney and Tremaine, 1987; 
Rahkov, 2003; LDT12). We handle the stirring and damping parts of dy¬ 
namical friction separately as shown below. 

The evolution of eccentricity e* and inclination ij of the target i due 
to viscous stirring and dynamical friction caused by interactions with the 
interloper j is given as (Ida, 1990; Tanaka and Ida, 1996) 


dej 

dt 

del 

dt 

del 

dt 


dil 


vs,i 


DFSj 


DFDj 


dt 


— C'l-PvS; 

'Hj 


C'l 


^jhlj 


2 e^Oj COSZUij^ Tj)F5 


- d^iQvs, 


VS,i 


DFSj 


'Hj 


( 17 ) 

(18) 

(19) 

( 20 ) 
( 21 ) 


dt 


DFD,j 


C'l 




Pdf ■) 


( 22 ) 


where Ci = njiyjhfjaljUJK, = nij/{mi+mj), Pys, Qys, and Pdf are the non- 
dimensional rates of viscous stirring and dynamical friction. The expressions 
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of these rates are given in Appendix B.2. Since we do not assume uniformities 
of Wij = Wj — Wi and ffjj = VLj — ffj to handle eccentric and/or inclined 
ringlets, as discussed in Section 2.3, the terms with these relative angles 
remain in the expressions for dynamical friction. The way to split dynamical 
friction into the two terms is explained in detail in Appendix C. 

The evolutions of due to viscous stirring and dynamical friction are 
then given as 



Similar expressions can be obtained for ii. 

Mutual encounters between planetesimals causes random walk in their 
semimajor axes. The radial diffusion coefficient Di for the target i is approx¬ 
imately given by (Appendix D) 


Di — ^ \Cial,j{PYs + Qvs)] • (25) 

j 


Ignoring the effect of curvature, the probability function of the change in a* 
after the time step At is given by 


P{Aai) 


1 ( 

V4irB.Ai ^ V -iDiAtJ' 


(26) 


where the function is normalized as P{Aai)dAai = 1. We choose Aa* at 
random with a weight given by P(Aaj). If Di < 0, we set Aa* = 0. This is 
rare but can occur if 3> fr- 


2.5.2. Conversion of stirring rates to accelerations 

The stirring and damping rates of orbital elements need to be converted 
into the acceleration of the tracer. Since the time step for the statistical 
routine At is usually much smaller than the stirring time, ef {dej/dt)~^, the 
absolute values of the changes in the orbital elements during At are much 
smaller than the magnitudes of the orbital elements themselves. Consider a 
situation that encounters between planetesimals occur randomly in directions 
of velocities. A certain acceleration occurs to a planetesimal in an encounter. 
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and in another encounter an acceleration with the same magnitude to the first 
encounter but in the opposite direction can occur. In first order, the changes 
of orbital elements cancel out in these two encounters. However, in second 
order, the orbital eccentricity and inclination increase a little on average, be¬ 
cause the orbital elements are slightly changed after the hrst encounter. This 
is the way how viscous stirring works, and to handle it, we introduce second 
order Gauss planetary equations. The stirring part of dynamical friction is 
also handled by these equations. On the other hand, the damping part of 
dynamical friction is handled by the standard Gauss planetary equations of 
hrst order (e.g., Murray and Dermott, 1999). If we use the standard Gauss 
planetary equations for the stirring part of dynamical friction, the standard 
deviations of orbital eccentricity and inclination become unphysically small 
for equal-mass planetesimals or mass-dependent segregation of these orbital 
elements occurs for a case of a mass distribution. We assume that e* 1 
and 1 in the following formulation. Note that these assumptions are 
already used in Hill’s approximation. 

In this subsection only, we drop the subscript i for the tracer i to avoid 
long, confusing subscripts. The expected changes of orbital elements dur¬ 
ing the time step At are Aa, Acys + Ae^ps, Acdfd, Aiyg -|- Ai^pg, and 
AiDFD, where Acyg = {de^/dt)YsAt, Aeppg = {de^/dt)T)FsAt, Acdfd = 
(e^ -|- ((ie^/(it)DFDAt)^/^ — e, and similar expressions are obtained for the 
inclination. 

The change of the velocity vector is split into three components as 

Av = Av-rXji + Aut^t + Aun^^n, (27) 

where Aur and Aux are the components of the velocity change in the radial 
and tangential directions in the orbital plane, Aun is the velocity change 
normal to the orbital plane, a:;R, xt, and ojn are unit vectors. The velocity 
vector of the tracer itself is given as v = URasR -|- vtXt, where ur and ux are 
the radial and tangential components. All the components of the velocity 
change are further split into the stirring and damping terms with subscripts 
S and D, respectively: 

Aur = Aur^ -I- Aur^d, 

Aux = Anx,s T Anx,D) (28) 

Aun = AnN,s + Aun.d- 

In our stirring routine, all the components of the stirring part are chosen 
using random numbers, but their averaged functional forms follow Gaussian 
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distributions with the zero mean values: (Aur^s) = (Aut,s) = ('^'?^n,s) = 0, 
where the angle brackets in this subsection mean averaging over multiple 
velocity changes. The acceleration due to an individual encounter may not 
necessarily follow Gaussian. On the other hand, since we handle the time 
averaged accelerations of tracers due to multiple encounters that occur in a 
random fashion, their distributions are limited to Gaussian. 

First we consider the velocity changes in the orbital plane. The energy 
conservation gives the following equation, if we take up to the second order 
terms 

fi { Aa ( Aa\^\ 
a \ a \ a J I 

2urAur + 2utAut + (Aur)^ + (Aut)^ + (Aun)^, (29) 


where p = If we ignore the second order terms and the hrst term in 

the r.h.s. of Eq. fl2^ as ur -C ut, this gives Aut,s as 


Aut,s = 


pAa 

2a‘^VT' 


(30) 


Note again that Aa is given by Eq. fl26|l and thus Aut,s is Gaussian. The 
angular momentum conservation gives 


—a/i(2eAe + (Ae)^) + Aa/i(l — — 2eAe) = 

(2utAut + (Aut)^ + (Aun)^) , (31) 


where Ae is the change of e and R is the distance from the central star. 
Inserting Eq. fl^ into Eq. fl2^ . we have 


/iAa 


^ Aa a? o ^ A N 


+ 


fia{2eAe + (Ae) 


2urAur + (Aur)^ 


(32) 


Averaging Eq. fl5^ over 
only the second order terms. 


time, the hrst order terms vanish and we have 
Using Ae^ = 2eAe + (Ae)^, we have 


^(Ae2) = + ((Aur)2) 

a O'* 

= —(2^evs + 2^^vs) + ((^^^'r)^)) 
a 


( 33 ) 
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where we used a? ~ and ((Aa)^)/2 ~ a^(Aeyg + Aiyg) [~ -DjAt; Eq. fl2B]) 
or Eq. fl99p ]. and ignored the term proportional to AaAe as it is small enough 
as compared with other terms for e 1. Inserting (Ae^) = Acyg + Ae^pg 
into Eq. ([33]), we have 

((Aur^s)^) = ^ (Ae^pg — Acyg — 2Aiyg) . (34) 

(Jj 

We give Aur^s assuming it is Gaussian. If the r.h.s side of Eq. (IMll is negative, 
we set Aur^s = 0, add the r.h.s divided by p/a to /dt)BFD^t, and re¬ 
calculate Acdfd; which will be handled as shown below. This occurs in the 
shear dominant regime [(e^ -f < 2hjj]. 

To give the velocity changes Aur r and AnT,D for the damping parts of 
dynamical friction, we use the standard first order Gauss planetary equation 
for e (Murray and Dermott, 1999) that can also be derived from Eq. fl32p : 


Ae 




[AUr sin / -|- AnT(cos / -|- cos E )], 


(35) 


where / and E are the true and eccentric anomalies. We employ the following 
forms ^ 

Aur^d = Auro sin /, Aut.d = cos /. (36) 

Setting Ae = Acrfd in Eq- (|3S|), we have 


^^Ro — . , ^ —^Acdfd, (37) 

Y a(l — e^) 

where we assumed cos / ~ cosE. 

Next, we consider the velocity change perpendicular to the orbital plane. 
The specific angular momentum vector h = {h^, hy, h^) after the velocity 
change Aun in the inertial Gartesian coordinate system are given as 

(K\ ( 7?sin/AnN \ 

hy = P3P2-P1 -Pcos/Aun , ( 38 ) 

\h.) \{h^ - 

where h = \h\, P 3 , P 2 , and Pi are the rotation matrices given by Eqs. (2.119) 
and (2.120) of Murray and Dermott (1999). Using hz = hcos {i + Ai), the 
z-component of Eq. (I38P is given as 

Azsm*-|--—cos* =-—Aun cos (tc-|-/) smz-|---^(A un) cos*, (39) 

^ ! L ^ 11 
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where Ai is the change of i and w = w — is the argument of pericenter. 
This equation leads to the standard Gauss planetary equation for i (Murray 
and Dermott, 1999), if we ignore the second order terms. Averaging Eq. fl3^ 
over time, the hrst term in the r.h.s. vanishes. Setting the mean change of 
as (Ai^) = Aiyg + Aippg and using i -C 1, we have 

{(AnN,s)^) = + "^^DFs)- (40) 

CL 

We give Awn.s, assuming it is Gaussian. For the damping part of dynamical 
friction, we adopt the form 

Aun.d = Ar;NoCos(w +/). (41) 

Ignoring the second order terms in Eq. ([39]), and setting Ai = Aiopo, we 
have 

A 2h. ^ ^ , 

Avno = —Azdfd, (42) 

a 

where we took averaging over the orbital period and used J cos^ {w + f)df / (27r) 
1 / 2 . 

2.6. Treatments of sub-embryos 

If the planetesimal mass is equal to or exceeds mto, the tracer is 
promoted to a sub-embryo. The acceleration of a sub-embryo due to gravita¬ 
tional interactions with a single planetesimal in a tracer is taken into account 
in the A-body routine to hold the correct mutual Hill radius. On the other 
hand, the acceleration of the sub-embryo due to interactions with other kj — 1 
planetesimals in the tracer is handled by the statistical routine. Gollisions of 
planetesimals with sub-embryos are also handled in the statistical routine. A 
large embryo tends to form a gap around in its feeding zone (Tanaka and Ida, 
1997), although other nearby embryos, if they exit, tend to £11 the gap. The 
surface density in the feeding zone usually decreases through this process. 
We take into account a possible difference of surface densities inside and out¬ 
side of the feeding zone whereas spatial variation of the surface density inside 
of the feeding zone is not considered. This assumption allows us to still use 
the stirring and collision rates averaged over phases and semi-major axis for 
interactions between a sub-embryo and tracers. 
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2.6.1. Collisions and gravitational stirring 

A collision between the snb-embryo i and the tracer j occnrs only if 
E] > 0, where Ej is the Jacobi integral of Hill’s eqnations given by 




g {{dEoy - {dijf) Wk, 


(43) 


where dij = Oj — a* and is the half width of the feeding zone given by 


deo — 


V3 K, 



1/2 


(44) 


Even if Ej > 0, collisions do not occnr with trojans or if \dij\ is too small 
(see Fig. 2 of Ohtsnki and Ida, 1998). Ida and Nakazawa (1989) showed that 
collisions occnr only if 


\dij\ > dmin = [8/(2.5 + 2.0eJhij)Y^^aihij, (45) 


where dmin is the minimum orbital separation for occurrence of collision. We 
do not consider interactions with trojans in the statistical routine. Planetes- 
imals in the feeding zone {Ej > 0 and \dij\ > dmin) can potentially enter the 
Hill sphere of the embryo and orbits of planetesimals are strongly deflected 
by close encounters with the embryo. 

For tracer-tracer interactions, the number density is calculated by Eqs. o 
and (E]). If a sub-embryo is the target, however, we search for planetesimals 
in the feeding zone. Their surface number density is given by 

h ■ T- 

n- = _ ^ _ -2- 1461 

' 4r,(dso-c?min)ddT,-i,- ^ ^ 

We miss some of tracers in the feeding zone if they are radially outside of the 
region i of neighboring search. The factor Tj/Tj^in in Eq. fj46l) is introduced to 
correct this effect, where Tj is the orbital period of the tracer j, and Tj^in is the 
period in which the tracer j is radially in the region i during Tj provided that 
the tracer j is azimuthally in the range of the region i. The expression of in 
is given in Appendix F. In the stirring routine, we use kj — 1 in Eq. fl46l) instead 
of kj as the back reaction from a single planetesimal is taken into account in 
the Wbody routine. Using Eq. (061), collisional growth and velocity changes 
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of sub-embryos are handled in the statistical routine in the same manner as 
tracer-tracer interactions. If there is no spatial density variation between the 
inside and outside of the feeding zone, the same collision and stirring rates 
are derived whether the target is a tracer or a sub-embryo. 

Even if the mutual distance between the centers of a tracer and a sub¬ 
embryo is less than the sum of the radii of the constituent planetesimal and 
the sub-embryo, orbital integrations of both bodies are continued in the N- 
body routine using softened forces that avoid singularity. We employ the K1 
kernel of Dehnen (2001) for the softening. 

Gravitational interactions occur even if Ej < 0 due to distant encoun¬ 
ters (Weidenschilling 1989, Hasegawa and Nakazawa, 1990). We evaluate 
the change rates of the orbital elements of the sub-embryo due to individ¬ 
ual distant encounters without assuming uniformity in dij. Let the changes 
of the squared relative orbital eccentricity and inclination due to a single 
distant encounter be 6e^ and Hasegawa and Nakazawa (1990) derived 
the changes {Se^) and {6i‘^) averaged over and Hr [their Eqs. (38) and 
(39)]. The averaged changes of ej and of the sub-embryo i during a single 
encounter with a planetesimal in the interloper j are given by z/|(5e^) and 
(Ida, 1990). Using the Keplerian shear velocity given by (Bromley 
and Kenyon, 2011) 



(47) 


the frequency of encounters is given as VshnjSr = Vsh{kj — l)/{4:ai69), where 
we used Eq. ([2]) and 6r in this case is the radial extent of planetesimals in 
the tracer j and is assumed to be inhnitesimally small. The change rates of 
and due to distant encounters with the interloper j are then given as 



(48) 


(49) 


where Ri = 0.747 and = 0.210. The rates given by Eqs. flT8|) and fl49l) 
are added to Eqs. ffT7|l and fl20l) if Ej <0. Dynamical friction due to distant 
encounters is negligible (Ida, 1990). 
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2.6.2. Planetesimal-driven migration 

An embryo gravitationally scatters surrounding planetesimals. Back re¬ 
action on the embryo causes its radial migration. Individual gravitational 
encounters usually cause a random walk of the embryo. However, if the ac¬ 
cumulated torque is not cancelled out, the embryo migrates over a long dis¬ 
tance. This is called planetesimal-driven migration (Ida et ah, 2000; Kirsh 
et ah, 2009). Planetesimal-driven migration of full-embryos is automatically 
taken into account in the iV-body routine. We developed a unique method 
to handle planetesimal-driven migration of sub-embryos explained in the fol¬ 
lowing. 

In our code, the angular momentum change of the sub-embryo 

i due to gravitational interactions with tracers during At is stored in the 
iV-body routine. If AL^^i is released instantaneously by giving the tangen¬ 
tial acceleration, the sub-embryo may migrate unrealistically too fast due to 
strong kicks from tracers that have masses similar to the sub-embryo. The 
direction of migration may change rapidly in such a situation. This prob¬ 
lem can be avoided by limiting the angular momentum released during At 
to the theoretical prediction, ALs,i, derived from the statistical routine. If 
\AL]sf^i\ < \ALs,i\, AL^^i is fully released. If \ALN^i\ > lAL^^il, on the other 
hand, the remaining portion ALjq^i — ALs^i is added to AL^r^j stored during 
the next time step interval. What is done by the scheme is basically time 
averaging of the torque on the sub-embryo over a timescale much longer than 
At so that strong spikes of the torque due to close encounters are smoothed 
out. More explanations are given in Section 3.3, together with test results. 

The predicted angular momentum change, ALs,i, is determined from the 
surface density and velocities of planetesimals in neighboring tracers. Brom¬ 
ley and Kenyon (2011, 2013) and Ormel et ah (2012) showed three different 
modes of planetesimal-driven migration, referred to as Type I - III migration 
analogous to gas-driven migration. Type I and III migration occurs when an 
embryo interacts with planetesimals in close encounters. Type I migration 
is considered to occur when the velocity dispersion of planetesimals is large, 
and its migration rate is very slow as is determined by the difference between 
torques from outer and inner disks. On the other hand. Type III migration is 
fast self-sustained migration that occurs when the velocity dispersion of plan¬ 
etesimals is relatively small, and the embryo migrates due to the torque from 
the one side (inner or outer side) of the disk, leaving a relatively enhanced 
disk in the other side (Ida et ah, 2000; Kirsh et ah, 2009; Levison et ah, 2010; 
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Capobianco et al., 2011). We ignore Type I migration in the evalnation be¬ 
low, as the migration rate is well represented by the formnla for Type III 
migration even at large velocity dispersions (Kirsh et al., 2009; Bromley and 
Kenyon, 2011). Type II migration occnrs when an embryo opens np a gap in 
the planetesimal disk so the torqne is primarily from planetesimals in distant 
enconnters. This torqne is particnlarly important if the asymmetries of the 
snrface density and the velocity dispersion between the inner and onter disks 
are large. Snch a sitnation occnrs if an embryo in a gaseous disk shepherds 
outer small planetesimals that rapidly spiral toward the central star without 
the embryo. 

The radial migration rate of Type III planetesimal-driven migration due 
to an interaction with the tracer j is (Ida et ah, 2000) 



(50) 


where the sign is positive if aj > ai and vice versa, and Rug is the non- 
dimensional migration rate given in Appendix B.3. If ALjsf^i stored in the 
Wbody routine is positive/negative, we count only tracers with semimajor 
axes larger/smaller than a* because planetesimals in close encounters cause 
attractive forces. In the evaluation of rij in Eq. fl50|) . we use 0.5S'i(l±0.55r/rj) 
instead of Si used in Eq. (|2]). The Type III migration rate is then given as 



(51) 


where /slow is the factor that represents slowing down for massive embryos. 
The slowing down occurs when the embryo cannot migrate over the feeding 
zone during the synodic period for that width, as cancellation of angular 
momentum exchange occurs in the second encounter (Ida et ah, 2000, Kirsh 
et ah, 2009). Bromley and Kenyon (2011) showed that the slowing down 
factor is given by a migration length during the synodic period for the half¬ 
width of the feeding zone d{ relative to df, if this ratio is less than unity: 



(52) 


We employ df = l.Sajhy after Bromley and Kenyon (2011). This length 
is somewhat shorter than the actual half-length of the feeding zone, (Ieo 
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[Eq. dS])], and the reason for this valne being applied is becanse the strong 
attractive force conies from a relatively narrow region. If the embryo mass 
is small enongh, /glow = 1; in this regime, the migration rate is independent 
of the embryo mass. 

One may think that the migration rate can be simply estimated by snm- 
ming np the torqnes from both sides in Eqs. flSU]) and fl3T]l rather than snm- 
ming np for one side and that AL^v,* calcnlated in the A-body rontine is 
nnnecessary. However, we hnd that snch an approach can work only if the 
embryo is massive enongh so that /glow <1- In this case, the all tracers in 
the feeding zone strongly interact with the embryo and this prodnces large 
velocity asymmetry between the nnpertnrbed side and the other side that 
the embryo migrated throngh. On the other hand, if the embryo is small, it 
migrates over the feeding zone throngh interactions with a small fraction of 
tracers in the feeding zone, and only small velocity asymmetry is prodnced. 
In snch a sitnation, the evalnated torqnes from both sides ronghly cancel ont 
each other and the fast migration cannot be reprodnced. If we evalnate the 
torqnes from tracers in an azimnthally limited region, where velocity asym¬ 
metry exists, snmmation of both sides may work. However, an appropriate 
choice of the azimnthal width is nnclear to ns and the nnmber of tracers in a 
narrow region is too few to give a reliable migration rate. This is why we eval¬ 
nate the migration rate dne to the one-sided torqne as an npper limit, while 
the actnal angnlar momentnm change is evalnated in the A-body rontine. 

For Type H migration, we evalnate the torqnes from tracers in both sides 
of the disk as done for the stirring rate dne to distant enconnters. The 
change in the semimajor axis difference, 6dij, dnring a single enconnter for 
Cr = h = 0 is given as (Bromley and Kenyon, 2011) 



(53) 


where g{b) = —32b~^ and b = dij/{aihij). This formnla is expected to be inde¬ 
pendent of Cr and h after averaged over w,- and Hr (Hasegawa and Nakazawa, 
1990). In a similar manner of Eq. (l48ll . the change rate of the semimajor axis 
of the embryo is given as 
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( 54 ) 



The Type II migration rate is then given as 



( 56 ) 


If the surface density is symmetric for the inner and outer disk, the first term 
in the parentheses in the r.h.s. of Eq. (I5T|) vanishes after summation. In this 
case, the embryo slowly migrate inward. If density asymmetry exists, the 
hrst term remains and relatively fast migration takes place in either inward 
or outward direction. 

The migration rate of the embryo suggested from the statistical routine 
is given as 



(56) 


The predicted angular momentum change during At is given as 



(57) 


where L, is the angular momentum of the sub-embryo. To produce smooth 
Type III migration, at least a few, preferably ~ 10, tracers need to exist 
constantly in the feeding zone in the region i. If we choose a very small 66 
for the region z, that makes sub-embryo migration noisy. 

We hnd that a large AL^^i sometimes remains unreleased even when 
an embryo reaches a disk edge. This usually causes the embryo migrate 
artihcially further away from the disk edge. To avoid this problem, we remove 
quickly on the timescale of 1000 orbital periods at the disk edge, 
although this violates conservation of the total angular momentum of the 
system. We regard the embryo being at the disk edge if the total mass of 
tracers in the feeding zone ahead of the embryo in the migrating direction is 
less than the embryo mass. For a small embryo, sometimes no tracer exists 
in the feeding zone due to statistical fluctuations even if the embryo is not 
at the disk edge. To avoid this issue, we also count tracers within lOajhjj 
ahead of the embryo even if they are not in the feeding zone. 

2. 7. Global gravitational force of tracers 

The global gravitational forces of tracers are added to the equations of 
motion for tracers and sub-embryos. To calculate the forces, we assume that 
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the disk is axisymmetric and the vertical distribution is given by Gaussian 
for tracers. For simplicity, at each radial location, we consider only one scale 
height of tracers with various planetesimal masses whereas in reality different 
mass groups have different scale heights. The scale height of planetesimals 
in a radial grid at r is given by 



where the summation is done for each radial bin, and the total mass of the 
tracers Mtr in the grid is 


Mtr{r) = 'Yh 

tracers 


(59) 


The surface density of tracers are given by Str(r) = Mtr(r)/( 27 rrAr), where 
Ar is the radial grid width. The spatial density of tracers is given as 

The r and components of the tracers’ global gravitational force are given 
by (Nagasawa et ah, 2000) 
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where G is the gravitational constant, K{A and E{A are the first and second 
kind of elliptic integrals, and ^ is given by 
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These components of the force are tabulated in (r, z) grids and interpolated 
for use in simulations. For the accelerations on sub-embryos, — 1 is used 
instead of ki in Eqs. fl58|l and fl5^ . It is not necessary to calculate the global 
gravitational force components frequently and we usually update them every 
~ 1000 times the shortest orbital period of particles. 

The global gravitational force modihes the precession rates of the longi¬ 
tudes of pericenter and ascending node; it usually accelerates the absolute 
precession rates. The precession rates do not usually matter to dynamics of 
tracers in terms of such as encounter velocities. These rates are, however, 
likely to be important for stability of the system, when multiple sub-embryos 
appear in the system. The embryos mutually interact through secular reso¬ 
nances which cause oscillations in eccentricities and inclinations of embryos. 
The global gravitational force reduces these oscillation periods. As a result, 
the maximum eccentricities and inclinations of sub-embryos during the os¬ 
cillations are suppressed to relatively small values. We identify this effect in 
simulations in Section 3.5 (results without the global gravitational force are 
not shown there, though). Although its effect is much smaller than dynamical 
friction, it does not seem to be negligible. 

3. Tests of the hybrid code 

We show six different types of tests of the hybrid code: collisional damp¬ 
ing and stirring (Section 3.1), gravitational scattering and radial diffusion 
(Section 3.2), planetesimal-driven migration of sub-embryos (Section 3.3), 
runaway to oligarchic growth (Section 3.4), accretion of terrestrial planets 
from narrow annuli (Section 3.5), and accretion of cores of jovian planets 
(Section 3.6). The parameters used in the test simulations are summarized 
in Table 1. 

3.1. Collisional damping and stirring 

The first test we present is velocity evolution of tracers due to mutual 
collisions only to check whether the collision probability and the velocity 
changes at collisions are calculated accurately. All mutual-gravity related 
effects are ignored in the test: gravitational scattering, gravitational focusing 
in the collisional probability [the second term in the bracket of the r.h.s of 
Eq. fl69|l ]. and enhancement of the impact velocity (Appendix E.2). Collisions 
are assumed to result in inelastic re-bounds (hit-and-run collisions), and post 
collision velocities are characterized by the normal and tangential restitution 
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Section 

Theme 

Figure 

3.1 

Collision 

5 

3.2 

Scattering 

6-8 

3.3 

Migration 

9,10 

3.4 

Runaway 

11,12 

3.5 

Terrestrial 

13,14 

3.6 

Jovian 

15,16 

A^Tot (M®) 

0.1 

0.1,0.125 

230 

0.67 

2 

200 

Region (AU) 

0.95-1.05 

0.95-1.05 

14.5-35.5 

0.915-1.085 

0.7-1.0 

4-16 

a 

-1.5 

-1.5 

-1.0 

-1.5 

-1.5 

-1.5 

N 

o 

1—1 

1700-2 X 10® 

2.3 X 10® 

8.33 X 10® 

2000 

lOi® 

Ntr 

500,2000 

200,250 

230 

400,4000 

200,400 

10000 

Nem 

0 

0 

1 

> 0 

> 0 

> 0 

(e2)i/2 

0.005,0.0005 

io-‘‘,io-® 

0.01-0.05 

1.4 X 10-4 

0.001 

10-4 

St (days) 

6 

6 

60 

6 

6 

60 

SO (tt) 

0.5,0.1 

0.5 

0.5 

0.5,0.25 

0.5 

0.25 

Collision 

Rebound 

N/A 

N/A 

Merging 

Merging 

Merging 

Gas disk 

No 

No 

No 

Yes 

No 

Yes 


Table 1. Parameters used for hybrid simulations. Theme is a very brief ex¬ 
pression of what experiments are performed. Mpot is the initial total mass 
of planetesimals. Region represents the initial inner and outer edges of the 
planetesimal disk, a is the power-law exponent for the surface density of 
planetesimals (oc r“). N is the initial number of planetesimals. is the 
initial total number of tracers. iVem is the initial number of embryos (> 0 
means that no embryo is introduced in the initial state but embryos appear 
as a result of collisional growth), is the initial root-mean-square ec¬ 

centricity of tracers. St is the time step for orbital integrations. S6 is the 
angular half-width of the region used in neighboring search (Fig. 2). Collision 
represents collisional outcome. Gas disk represents whether damping effects 
due to a gas disk are taken into account. In all the simulations, the physical 
density of any body is 2.0 g cm“^. The initial distributions of e and i are 
given by the Rayleigh distributions with 
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0.005 



(b) 



Fig. 5. Time evolution of and of planetesimals due to mutual 

collisions only for (a) en = 0 and (b) = 0.8. All effects related to mutual 
gravity are ignored. The total number of planetesimals is 10^^ represented by 
2000 (solid curves) and 500 (dotted curves) tracers. Planetesimals are placed 
between 0.95 AU and 1.05 AU and the total mass is 0.1 M®. In the hybrid 
simulations, we only use tracers between 0.96 and 1.04 AU in calculations 
of and to avoid the effect due to radial expansion. Evolution 

curves calculated using the analytic formula (Ohtsuki, 1999) are also shown 
by dashed curves for comparison. In the panel (b), Cn = 0.72 is used for the 
analytic curves. 
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coefficients, Cn and Ct- Appendix E.2 explains how to handle hit-and-run 
collisions between planetesimals in two tracers, with assumptions of Cn = 0 
and Ct = 1 [see Eq. fll02p for more detail]. These coefficients, however, 
can be modihed to any values. In the tests, we assume no friction (i.e., 
Ct = 1). Collisional velocity evolution is a well-studied problem for planetary 
rings and several authors have developed analytical formulae (Goldreich and 
Tremaine, 1978; Hameen-Antilla, 1978; Ohtsuki, 1992, 1999). We compare 
our simulation results with analytic estimates of Ohtsuki (1999) because of 
the easiness of use. 

Figure 5 shows the time evolution of and of 10^^ planetes¬ 

imals for (a) en = 0 and (b) Cn = 0.8. The radius of all planetesimals is 0.4 
km. The planetesimal disk extends from 0.95 to 1.05 AU and its total mass 
is 0.1 M^. These parameters give the disk normal optical depth of 5 x 10“®. 
Two simulations are performed to see resolution effects for both Fig. 5a and 
5b: one with Wr = 2000 and 59 = O.lvr and another one with Wr = 500 
and 59 = 0.57r, where Wr is the total number of tracers. The total angular 
momentum, is well conserved for all test simulations shown in Fig. 5; 

l-htot(^) ~ Liotit = 0)|/Ltot(^) < 10 

For Cn = 0 (Fig. 5a), we hud excellent agreements between our results 
and the analytic estimates of Ohtsuki (1999), in particular for the case of 
Atr = 2000. However, the decreases in and stall at certain 

values in hybrid simulations - ^ 10“‘^(4 x 10“'^) for Atr = 2000 (500) - 

whereas in the analytic curves, e and i decrease down to ~ s/r, where s is the 
particle radius. In the hybrid simulations, eccentric ringlets form, in which 
the longitudes of pericenter of all particles are nearly aligned. This issue is 
probably related to the initial condition. In creation of the initial positions 
and velocities, the longitudes of pericenter are randomly chosen, and the sum 
of the eccentricity vectors of all tracers is not exactly zero. The eccentricities 
of the ringlets in the end state are, in fact, close to the expected values of 
the residual eccentricity, = 0)/A|-y^. Since the residual eccentricity 

decreases with increasing the number of tracers, the eccentricity of the ringlet 
in the end state decreases. Even if the residual eccentricity is hxed to be zero, 
the same problem is likely to occur unless all particles effectively interact with 
each other. The ringlets are not only eccentric but also inclined because a 
similar discussion is applied to inclination. The evolution curves of e and 
i also weakly depends on 50; a small 59 gives a better agreement with the 
analytic curves. It is found that relative velocities are slightly underestimated 
if 59 is large and radial excursions of particles are larger than 5r. This leads 
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a slightly lower damping rate. 

Collisions cause not only damping but also stirring as the collisional vis¬ 
cosity converts the tidal potential to the kinetic energy. If €n > Cent = 0.627, 
where Cent is the critical restitution coefficient, the stirring rate exceeds the 
damping rate and and increase with collisions (Goldreich and 

Tremaine, 1978; Ohtsuki, 1999). We hnd that our method is slightly more 
dissipative than the analytic estimate and gives Ccrit — 0.7 for reasons that 
we were not able to hgure out. Figure 5b shows the same as Fig. 5a but for 
the case of = 0.8. It is found that our result for = 0.8 and Wr = 2000 
(500) well agrees with the analytic estimate with = 0.72 (0.70). Overall, 
our method can reproduce analytic results using a ~ 10% higher Cn. In sim¬ 
ulations of planetary accretion, we use Cn = 0 for hit-and-run collisions and 
the offset in is not a problem. Note that orbital alignments do not occur 
for the case of Fig. 5b, because collisional stirring randomizes orbital phases. 

3.2. Gravitational scattering and radial diffusion 

The next test we present is velocity evolution of tracers due to mutual 
gravity only in the absence of any collisional effects. Figure 6 shows evolution 
of and for a planetesimal disk extending from 0.95 and 1.05 

AU. The panels (a) and (b) are for mono-sized planetesimals represented 
by 200 tracers and the numbers of planetesimals N are 2000 and 2 x 10®, 
respectively. The panel (c) is for a bimodal size distribution case; 1500 small 
planetesimals and 200 large planetesimals are placed in the disk. The mass 
ratio of a large planetesimal to a small one is 5. The number of tracers 
for small planetesimals is 150 and that for large planetesimals is 100. The 
results are compared with analytic estimates of Ohtsuki et al. (2002). For 
the panels (a) and (c), we also performed pure iV-body simulations using the 
parallel-tree code pkdgrav2 (Morishima et al., 2010) with the same initial 
conditions. The agreements between the analytic estimates and pure iV-body 
simulation results are excellent, so results from both methods are trustable 
as benchmarks. 

Overall, we hnd good agreements between the results from the hybrid 
code and the analytic estimates, although our model gives slightly lower 
and than the analytic estimates, except Fig. 6b shows a good 

agreement. We also perform the same hybrid simulations with different values 
of Ntr and 59, and the absolute levels of and at a certain time 

are found to be almost independent of these parameters, unless Wr and 59 
are too small. The underestimates of and of our calculations 
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(a) 



(b) 



Fig. 6. 
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(c) 



Fig. 6. Time evolution of and of planetesimals due to gravita¬ 

tional interactions: (a) N = 2000, A^"tr = 200, M^ot = O.IMq, (b) N = 2x10®, 
Nt, = 200, Mtot = O.lMe, and (c) N = 1700, = 250, Mtot = 0.125Me 

(see Table 1 for the characters). The panel (c) is the case of the bimodal 
size distribution with Ni = 1500, N 2 = 200, iVti -1 = 150, and iVti. 2 = 100, 
where the characters with the subscripts 1 and 2 are those for small and large 
planetesimals, and the mass ratio of a large planetesimal to a small one is 
5. In all cases, planetesimals are placed between 0.95 AU and 1.05 AU, and 
and are calculated using tracers between 0.96 AU and 1.04 AU. 

Evolution curves calculated using the analytic formula (Ohtsuki et ah, 2002) 
and from the pure 77-body simulations (only for the panels (a) and (c)) are 
also shown for comparison. 


35 

















Eccentricity 


0.012 
0.010 
0.008 
0.006 
0.004 
0.002 
0.000 
0.012 

0.010 
0.008 
0.006 
0.004 
0.002 
0.000 
0.012 

0.010 
0.008 
0.006 
0.004 
0.002 
0.000 
0.012 
0.010 
0.008 
0.006 
0.004 
0.002 
0.000 

0.95 1.00 1.05 0.95 1.00 1.05 

Semi-mojor axis (AU) Semi-major oxis (AU) 



0.012 

0.010 

0.008 

0.006 

0.004 

0.002 

0.000 

0.012 

0.010 
0.008 
0.006 
0.004 
0.002 
0.000 
I 0.012 

uj 0.010 
0.008 
0.006 
0.004 
0.002 
0.000 
0.012 
0.010 
0.008 
0.006 
0.004 
0.002 
0.000 


160 yr 


Pure N-body 




800 yr 



4800 yr 





16000 yr 



Fig. 7. Snapshots on the a — e plane for the simulations shown in Fig. 6a 
(a) the hybrid simulation and (b) the pure iV-body simulation. 



Fig. 8. Time evolutions of the standard deviation of 1/a for the simulations 
shown in Fig. 7. Results for six individual hybrid simulations are shown by 
colored dotted curves, whereas their mean value is shown by a black solid 
curve. The result for the hybrid simulation shown in Fig. 7a is a black dotted 
curve. For comparison, the result from the pure N-body simulation (Fig. 7b) 
is shown by a black dashed curve. 

is probably caused by the underestimates of Pys and Qvs that we employ. 
The analytic calculations of Tanaka and Ida (1996) give slightly lower Pys 
and Qvs for e^/hij > 4 and ir/hij < 1 than those from direct three-body 
orbital integrations (Ida, 1990, Rafikov and Slepian, 2010). The influence of 
this difference is less important for large {e^y^'^/hij and and that 

is why the hybrid simulation in Fig. 6b shows an excellent agreement with 
the analytic result, contrary to Fig. 6a. 

Figure 7 shows snapshots on the a — e plane for the simulation of Fig. 6a. 
Snapshots of the corresponding pure A^-body simulation are also shown for 
comparison. The hybrid code can reproduce not only and but 

also the distributions of e and i. Ida and Makino (1992) carried out iV-body 
simulations and found that the distributions of e and i can be well approx¬ 
imated by a Rayleigh distribution. We perform the Kolmogorov-Smirnov 
(KS) test to check whether the distributions of e and i in the hybrid sim¬ 
ulations resemble a Rayleigh distribution. The test gives the signihcance 
level, Qks (0 < Qks ^ 1) (Press et ah, 1986). A large Qks rneans that 
the likelihood that two functions are the same is high. The time-averaged 
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(0 < t < 50, 000 yr) values of Qks for the distributions of e and i are 0.44 
and 0.57, respectively. We also perform the same KS test using the outputs 
of the pure iV-body simulation, and the values of Qks for e and i are 0.62 
and 0.19. Thus, the distributions of e and i from our hybrid simulations re¬ 
semble a Rayleigh distribution at a level similar to those from pure Wbody 
simulations. 

Figure 7 also shows that a planetesimal disk radially expands with time 
due to viscous diffusion at a rate similar to that for the pure iV-body sim¬ 
ulation. To compare more quantitatively, we plot the time evolution of the 
standard deviation of 1/a (the normalized orbital energy) of the disk in Fig. 8. 
Since the radial diffusion of a low resolution hybrid simulation is found to be 
noisy, we take average over several simulations. The off-sets between simu¬ 
lations at t = 0 are caused by random numbers that are used for choosing 
initial semimajor axes of tracers. We hnd a good agreement in the increasing 
rates of the standard deviations of 1/a between the hybrid simulations and 
the pure iV-body simulation. This comparison indicates that the diffusion 
coefficient given by Eq. fl25]) is a good approximation. 

Contrary to collisional velocity evolution shown in Section 3.1, conserva¬ 
tion of the angular momentum is violated by gravitational scattering as we 
do not explicitly solve pair-wise gravitational interactions in the statistical 
routine. Provided that the radial migration distance of a tracer due to ran¬ 
dom walk is about ajCj, the error in the total angular momentum is estimated 
as |Ltot(^) — -htot(^ = 0)|/Ltot(^) ~ ■, which we found consistent 

with the actual errors in the test simulations. 

3.3. Planetesimal-driven migration of sub-embryos 

In Section 2.6.2, we described how we handle migration of sub-embryos. 
To test the routine of sub-embryo migration, we perform a suite of hybrid 
simulations starting with initial conditions similar to those in Kirsh et ah 
(2009) and Bromley and Kenyon (2011). The initial planetesimal disk ex¬ 
tends from 14.5 AU to 35.5 AU and its total mass is 230 represented by 
230 tracers (thus rujo = A sub-embryo with a mass of is initially 

placed at 25 AU and growth of the embryo is turned off. In this subsection, 
the index i is given to a sub-embryo. Each tracer has 10000 planetesimals 
and tracer-tracer interactions are quite small compared with gravitational 
effects of the embryo. 

An example of migration of an Earth-mass embryo {rui = M^) is shown in 
Fig. 9 using four snapshots on the a-e plane. For this example, the migration 
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Semi-mojor oxis (AU) 


Fig. 9. Self-sustained (Type III) migration of a sub-embryo. We em¬ 
ploy rrii = rrito = IM 0 (where rrii is the embryo mass), Wr = 230, and 
= 2{i'^y^'^/hij = 1. The embryo is displayed as a red hlled circle 
with a horizontal bar with a half length of 10 Hill radii. 
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(°) 



(b) 



(c) 



Fig. 10. Time evolution of the semi-major axis of a sub-embryo due to self- 
sustained migration: (a) = 1 and ruj = IM 0 , (b) /hij = 5 

and mi = IM®, and (c) jh^ = 1 and = 2 OM 0 . For all cases, 

A^tr = 230 and ruto = IMb- Analytic migration rates (Bromley and Kenyon, 
2011) are shown by dashed lines. 
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is in the mass independent regime [/slow = 1; Eq. (IS2]) ]. The embryo hrst 
has some strong enconnters with tracers and a stirred region is prodnced. 
The embryo tends to migrate away from the stirred region dne to strong 
torques exerted by tracers in the unperturbed region. This process trigers 
fast self-sustained migration. During the migration, the embryo continuously 
encounters with unperturbed tracers on the leading side while leaving the 
stirred tracers in the trailing side. The fast migration continues until the 
embryo reaches the disk edge. 

Embryo migrations for three cases with different embryo masses and ini¬ 
tial velocity dispersions are shown in Fig. 10. In each case, six simulations 
are performed, using different random number seeds to generate initial or¬ 
bital elements of tracers. The migration rates suggested by Bromley and 
Kenyon (2011) are shown by dotted lines. They conhrmed good agreements 
between their formula and A^-body simulations. Overall, our method can 
reproduce migration rates suggested by Bromley and Kenyon (2011). How¬ 
ever, our method gives almost the same probabilities of inward and outward 
migrations if rrii ~ ruto, whereas A^-body simulations show that inward mi¬ 
gration is predominant (Kirsh et ah, 2009; Bromley and Kenyon, 2011). If 
the embryo is small, the direction of migration is determined by hrst few 
strong encounters, and once migration starts it is self-sustained regardless of 
its direction. In realistic planetary accretion, even if inward migration is pre¬ 
dominant, the embryo encounters with inner massive embryos and then turns 
its migration direction outward (Minton and Levison, 2014). Therefore, we 
believe that the artihcially high occurrence of outward migration is not a 
problem. For a large sub-embryo, inward migration is clearly predominant 
even in our hybrid simulations (Fig. 10c). 

To understand these trends and why our method works, we consider con¬ 
ditions necessary to produce smooth self-sustained migration. If planetes- 
imals are small enough compared with the embryo, smooth self-sustained 
migration of the embryo takes place (ordered migration). If planetesimals 
are large, on the other hand, the embryo suffers strong kicks and the di¬ 
rection of its migration may change rapidly in a random fashion (stochastic 
migration). Whether migration is ordered or stochastic is determined by the 
size distribution of planetesimals. This argument is similar to that discussed 
for planetary spins (Dones and Tremaine, 1993). 

The relative contributions of the ordered and stochastic components are 
estimated as follows. Consider that an embryo migrates over a distance Dai 
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through n encounters with surrounding planetesimals: 


Dai = n{iyjaihij6b)^ 


(63) 


where again Uj = rrij / {mi+rrij) and Sb is the change of b [= dij /(ajhjj)] during 
a single encounter. In the above, the angle brackets mean the averaging over 
encountered planetesimals. The expected value of {DaiY is given as 


{Dai)^ 




n{{h>jaihij5bY) + n{n 


nalh%{u^){{5bf) 


1 + 


1 ) (^ 7 ^ 2^27 ^ 6 ) , 

x^])) \xm) 


{Daj) 

dij 


, (64) 


where we assume that n S> 1 and that there is no correlation between Uj and 
6b to transform to the second row. If the hrst term in the r.h.s of Eq. fl6T|) is 
larger than the second term, the migration is stochastic rather than ordered. 
This means that the embryo is likely to migrate over Dai via random walk 
rather than self-sustained migration. If we assume that the ordered torque 
exerted on the embryo is the one-sided (Type III) torque, (Sb) is determined 
by the non-dimensional migration rate Rug (Appendices B.3). The value of 
{{6bX) is estimated using the non-dimensional viscous stirring rates, Pys cind 
Qys (Appendix D), and we consider contribution of planetesimals in one side 
assuming that large velocity asymmetry exists. Using these non-dimensional 
rates, we have {6b) / {{6bX) ~ Pmg/(-Pvs + Qys)- This ratio is ~ 0.12 for the 
shear dominant regime [(e^+ < 2hij] and ~ 0.5hij/er for the dispersion 

dominant regime (we assume Cr = 2U). Therefore, even with planetesimals 
as massive as the embryo {{vj)/{vf) — 2), the migration is likely to be self- 
sustained migration rather than random walk if the embryo migrates longer 
than ~ ai{Ahij -|- Cr). This length is comparable to the feeding zone width. 

In the above discussion, the planetesimal mass rrij can be replaced by 
the tracer mass kjrrij to discuss migration calculated in our A^-body routine. 
Therefore, our A^-body routine can reproduce self-sustained migration of the 
embryo even with tracers as massive as the embryo. This is true only if the 
ordered torque exerted on the embryo is as strong as the one-sided torque. 
This is not always the case. The torque becomes nearly one-sided if a strong 
velocity asymmetry of planetesimals exists around the embryo. If the embryo 
jumps over its feeding zone due to a strong kick from a tracer that is as 
massive as the embryo, the embryo cannot see velocity asymmetry at that 
radial location. Thus, the condition to apply the one-sided torque is that 
radial migration distance of the embryo due to a single encounter is much 
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smaller than its feeding zone. This condition generally turns out to be Vj <C 1. 
It is likely that the condition Vj < 1/100 or less recommended by various 
authors (Kirsh et ah 2009; LDT12; Minton and Levison, 2014) comes from 
this criterion. In our method, alternatively, instantaneous large changes of 
Oj of the sub-embryo are suppressed by limiting the migration distance of 
the sub-embryo to the theoretical prediction that is derived assuming that 
planetesimals are much less massive than the embryo. Therefore, our method 
can reproduce a long distance self-sustained migration. 

The condition to reproduce predominant inward migration is more severe 
because we need to consider migration of an embryo due to torques from both 
sides of the disk. In an unperturbed state, the total torque from both sides is 
order of hij of the one-sided torque in the shear dominant regime (Bromley 
and Kenyon, 2011) and is order of of the one-sided torque in the dispersion 
dominant regime (Ormel et ah, 2012). Thus, [hij + e^) is multiplied to the 
second term in the r.h.s. of Eq. fIMl) . Then, the condition that the migration 
is ordered when the embryo migrates over the feeding zone is roughly given 
by {hij + Cr) > The initial condition of Fig. 10c barely satishes 

this criterion provided that rrij is the tracer mass instead of the planetesimal 
mass. Our method indeed reproduces predominant inward migration. Once 
the embryo migrates over the feeding zone, velocity asymmetry is obvious 
and self-sustained migration immerses. 

3.4- Runaway to oligarchic growth 

Most of the authors who developed statistical planetary accretion codes 
(Weidenschilling et ah, 1997; Kenyon and Luu, 1998; Inaba et ah, 2001; 
Morbidelli et ah, 2009; Ormel et ah, 2010a; LDT12; Amaro-Seoane et ah, 
2014) tested their codes using Fig. 12 of Wetherill and Stewart (1993) as a 
benchmark. We also perform simulations starting with the same condition of 
Wetherill and Stewart (1993). Inaba et ah (2001) performed the same statis¬ 
tical simulations using their own coagulation code and the code of Wetherill 
and Stewart (1993) after modifying settings of these two codes as close as 
possible. Both results are shown in Fig. 10 of Inaba et ah (2001) and we 
compare them with our results. As the initial condition, 8.33 x 10® planetes¬ 
imals are placed between 0.915 AU and 1.085 AU. Individual planetesimals 
have an equal-mass of mo = 4.8 x 10^® g. The nebular gas is included and 
the gas/solid ratio is 44.4. See Table 1 of Wetherill and Stewart (1993) for 
more details. We take into account aerodynamic gas drag (Adachi et eh. 
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Fig. 11. Snapshots on the a — e plane for the hybrid simulations; (a) A^tr = 
400 and (b) Ntr = 4000. Tracers are displayed as colored circles and their 
colors represent the mass of a planetesimal as shown by the color bar. Sub¬ 
embryos are displayed as red hlled circles with a horizontal bar with a half 
length of 10 Hill radii [for the panel (b), only sub-embryos more massive than 
lOrritQ have a horizontal bar for consistency with the panel (a)]. A radius of 
a circle is proportional to (kirrii)^^^. 
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(o) (b) 



Fig. 12. Snapshots of the cumulative number and the horizontal velocity 
(= (5/8)i/2g jUk, where uk is the Keplerian velocity at lAU) for (a) Ntr = 400 
and (b) Ntr = 4000. The diagonal dashed line is the number of planetesimals 
represented by a single tracer with a mass of rrito for each run. The vertical 
dotted line represents rrito- For comparison, the results from (c) Wetherill 
and Stewart (1993) and (d) Inaba et ah (2001) are shown; both from Fig. 10 
of Inaba et al. (2001). 
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Fig. 12. Continued. 
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1976)@ and tidal damping (Papaloizou and Larwood, 2000), but not Type I 
gas-driven migration. The gas disk model is described in Morisliima et al. 
(2010) in detail. In our simulations, the planetesimal disk radially expands 
with time, while the solid surface density remains constant in the previous 
statistical simulations. However, the effect of radial expansion is small on 
the timescale of the simulation ( 10*^ yr). 

Figure 11 shows snapshots on the a — e plane for the two hybrid sim¬ 
ulations with different initial numbers of tracers, Ntr = 400 and 4000. In 
the early stage of evolution {t = 10,000 yr), large embryos grow very quickly 
and e is enhanced locally where large embryos exist. This stage is called the 
runaway growth stage. Once large embryos dominantly control the velocity 
of the entire planetesimal disk {t = 25,000 yr), growth of embryos slows 
down (Ida and Makino, 1993). This stage is called the oligarchic growth 
stage. Even though large embryos tend to open up gaps around their orbits, 
the gaps are not cleared out because the embryos migrate quickly due to oc¬ 
casional short-distance planetesimal-driven migrations and embryo-embryo 
interactions and because other embryos tend to hll gaps. During the oli¬ 
garchic growth stage, the orbital separation between neighboring embryos is 
typically 10 Hill radii as found by Kokubo and Ida (1998, 2000). Tanaka 
and Ida (1997) showed that gaps open up only if the orbital separation of 
neighboring embryos is larger than 20 Hill radii. Thus, formation of clear 
gaps is unlikely to occur in this stage. 

Figure 12 shows the mass and velocity distributions at different times 
for two hybrid simulations shown in Fig. 11. The mass spectra of the hy¬ 
brid simulations are smooth at the transition mass rrito at which a tracer 
is promoted to a sub-embryo (see also Appendix G). For comparison, the 
results from Wetherill and Stewart (1993) and Inaba et al. (2001) are also 
shown. The hybrid simulations agree well with these statistical simulations, 
although some small differences are seen. In the hybrid simulations, the 
number of planetesimals represented by tracers has a lower limit shown by 
a dashed line in Fig. 12. Thus, the tail of the mass spectrum in the massive 
side during the runaway growth is not correctly described. This causes a sys¬ 
tematic delay in growth of large bodies, as already discussed by LDT12. The 
delay is particularly clear in the low-resolution hybrid simulation at 10,000 


al. 


^We employ the gas drag coefficient of Dd = 0.5 instead of 2.0 used in Morishima et 
( 2010 ). 
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yr, although the delay is small. This resolution effect should be seen even 
in the high resolution hybrid simulation while we see a good agreement with 
the statistical simulations at 10,000 yr. This is probably because our stirring 
routine slightly underestimates eccentricities and inclinations (Section 3.2) 
and this effect cancels out with the resolution effect to some degree. Once 
the growth mode reaches oligarchic growth (t > 25,000 yr), the masses of the 
largest bodies in the low resolntion simulation almost catch np with those in 
the high resolntion simulation and those in the statistical simnlations. 

3.5. Accretion of terrestrial planets from narrow annuli 

The last section focused only on the early stages of planetary accretion. 
Fnrther, the damping effects dne to the gas disk make it difficult to judge 
whether the code works appropriately. Here we test our code for the entire 
stages of accretion withont the gas disk. We perform hybrid simnlations of 
accretion of terrestrial planets starting from a narrow annulns, similarly to 
those adopted in Morishima et al. (2008) and Hansen (2009), nntil comple¬ 
tion of accretion {t ~ 100 Myr). As the initial condition, 2000 eqnal-mass 
planetesimals (mo = lO^^M®) are placed between 0.7 AU to 1.0 AU, and the 
total disk mass is 2.OM0. The initial nnmber of tracers is 200. For compari¬ 
son, we also perform pnre A^-body simulations with the same initial condition 
nsing pkdgrav2. A hybrid simnlation takes abont three days using a single 
core whereas a pnre Wbody simnlation takes a few weeks using eight cores. 
As fonnd in Morishima et al. (2008), the evolntion is almost deterministic 
dne to efficient dynamical friction if the initial nnmber of planetesimals is 
large enongh (A^ > 1000). If the initial N is small while keeping the same 
total mass, the evolntion becomes rather stochastic and the final systems 
have relatively large diversity (Hansen, 2009). 

An example of the hybrid simnlations is shown in Fig. 13 together with 
an example of the pnre A^-body simnlations. Throngh the rnnaway growth 
stage, tens of embryos form [t = 0.1 Myr). With increasing masses of em¬ 
bryos, mntnal orbital crossings and impacts occur, and a few distinctively 
large embryos emerge {t= 1 Myr). In the late stage, the largest embryos 
sweep np small remnant planetesimals {t > 10 Myr). At the end of accre¬ 
tion, the largest three embryos have orbital eccentricities and inclinations as 
small as those for the present-day Earth and Venus due to dynamical friction 
of remnant bodies {t = 100 Myr). In the hybrid simulation, all the largest 
embryos are snb-embryos nntil the end of the simnlation. We perform three 
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Fig. 13. Snapshots on the a — e plane for the simulations starting from equal- 
mass 2000 planetesimals in a narrow annulus between 0.7 AU and 1.0 AU: 
(a) the hybrid simulation starting with 200 tracers and (b) the pure A^-body 
simulation. In the panel (a), tracers are displayed as colored circles, and 
the colors represent the mass of planetesimal. Sub-embryos are displayed as 
red hlled circles. Sub-embryos more massive than lOmio (= O-IM®) have a 
horizontal blue bar with a half length of 10 Hill radii. The same color code 
is used in the pure A^-body simulation. A radius of a circle is proportional 
to (for the pure A^-body simulation, ki = 1 for all particles). 
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Fig. 14. 
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(c) 



Fig. 14. Time evolutions of various quantities for three hybrid simulations 
(orange, brown, and red solid curves) and three pure A^-body simulations 
(blue, purple, and black dashed curves): (a) the mass-weighted root-mean- 
square eccentricity and inclination, ((me^)/(m))^/^ and ((mi^)/(m))^/^, for 
embryos and planetesimals, (b) the total nnmber of all bodies N, the total 
nnmber of embryos Nem, the total nnmber of tracers Ntr, (c) the total mass 
of embryos Mem,tot, fhe effective mass of the system (m^)/(m), and the mass 
of the largest body rrii. 
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hybrid simulations and three pure iV-body simulations in total, and all simu¬ 
lations show similar results, except one hybrid simulation produces only two 
distinctively large embryos. 

To see the accretion history more in detail, we plot various quantities as 
a function of time in Fig. 14. Overall agreements between the hybrid sim¬ 
ulations and the pure iV-body simulations are excellent. Good agreements 
between the hybrid and pure iV-body simulations suggest that our statistical 
routine works well even in the late stages of planetary accretion. Figure 14a 
shows the evolution of the mass-weighted root-mean-square eccentricity and 
inclination, Cm = ((me^)/(m))^/^ and im = for embryos and 

planetesimals. With increasing time, and v for both embryos and plan- 
etesimals increase primarily due to viscous stirring of largest bodies. After 
~ 10 Myr, Cm and im of embryos starts to decrease because small embryos 
with large e and i collide with large embryos with low e and i . Figure 14b 
shows the evolution of the total number of all bodies N and the total num¬ 
ber of embryos Agm- The number N^m has a peak around 1 Myr both in the 
hybrid and pure A-body simulations. The number of bodies which orbits are 
numerically integrated, Atr -|- Agm, increases by ~ 30% at maximum and the 
maximum value is taken around 1 Myr (Fig. 14b). 

Figure 14c shows the total mass of embryos, Mgm.tot, the effective mass 
of the system {w?)/(m), which is sensitive to masses of largest embryos, and 
the mass of the largest body, mi. It is found that the values of {w?)/ (m) and 
mi in the hybrid simulations are slightly lower than those in the pure A-body 
simulations at f ~ 1 Myr. During the giant impact stage roughly between 
t = 0.2 Myr and a few Myr, the contribution of embryo-embryo impacts 
to mass gain of embryos is signihcant, comparably to embryo-planetesimal 
impacts. The embryo-embryo collision rate sensitively depends on e and i of 
embryos. The hybrid simulations give slightly higher e and i (Fig. 14a) and 
thus lower growth rates of embryos than those in the pure A-body simulations 
during the giant impact stage. This is likely to be a resolution effect. In 
the low resolution hybrid simulations, the number of nearby tracers around 
embryos fluctuates and that makes dynamical friction somewhat noisy. We 
also perform the same simulations but starting with a larger number of tracers 
(Atr = 400) and hnd better agreements in (m^)/(m) with the pure A-body 
simulations. 
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Fig. 15. Snapshots on the a — e plane for the hybrid simnlation in the jovian 
planet region. The simnlation starts with eqnal-mass 10^^ planetesimals. 
Tracers are displayed as colored circles and colors represent the mass of a 
planetesimal. Sub-embryos (O.O 2 M 0 < m* < 2 M 0 ) are displayed as red 
hlled circles and full-embryos (m* > 2 M 0 ) are displayed as blue hlled circles. 
A radius of a circle is proportional to but enhanced by a factor of 

6 for embryos for good visibility. 
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Fig. 16. Time evolutions of (a) semimajor axes and (b) masses of embryos 
that are more massive than O.dM^ at 2 Myr. The same color is used for 
the same embryo in both panels. The embryos represented by blue, red, and 
purple curves are called Embryo 1, 2, and 3, respectively, in the main text. 
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3.6. Accretion of cores of jovian planets 

The final test we present is accretion of cores of jovian planets starting 
with a wide planetesimal disk. Levison et ah (2010) and Capobianco et ah 
(2011) showed importance of self-sustained outward migration of embryos, as 
they can grow rapidly via collisions with little-perturbed outer bodies. How¬ 
ever, these authors placed large embryos in the disk of small planetesimals 
in the initial conditions, and it is not clarified if outward migration plays a 
significant role in more realistic accretion scenarios. 

We set the initial condition similar to that employed by Levison et ah 
(2010), but without any embryos. We place 10^^ equal-mass planetesimals 
between 4 AU and 16 AU represented by 10000 tracers. The initial plan¬ 
etesimal mass mo is 1.19 x 10^^ g and the radius is 2.4 km. The total mass 
of planetesimals is 200 and this gives the solid surface density about 
seven times the minimum-mass solar nebula (Hayashi, 1981). The simula¬ 
tion is performed up to 2 Myr and it takes roughly five days to complete 
using a single core. The nebula gas effects are included in the same manner 
as the simulations in Section 3.4, and the gas-to-solid ratio is 56.7, which is 
the same as the minimum-mass solar nebula beyond the snow line (Hayashi, 
1981). We do not consider enhanced cross sections of embryos due to their at¬ 
mospheres (Inaba and Ikoma, 2003). The nebular gas is assumed to dissipate 
exponentially on the timescale of 2 Myr. Capobianco et al. (2011) showed 
that outward migration of an embryo is usually triggered in the nebular gas, 
if the planetesimal radius is in a certain range. The planetesimal radius we 
employ is in this range in the region inside of ~10 AU [see their Eq. (21)]. 
Note that, however, planetesimal radii increase with time in our simulation 
unlike those in Levison et al. (2010) and Capobianco et al. (2011). 

The snapshots on the a — e plane are shown in Fig. 15 and the time 
evolutions of the masses and the semimajor axes of six largest embryos at 
the end of the simulation are shown in Fig. 16. Several embryos appear 
near the inner edge of the disk after 0.1 Myr. After some embryo-embryo 
interactions, one of the embryos (Embryo 1) starts self-sustained outward 
migration. Since Embryo 1 encounters with unperturbed planetesimals on 
its way, it grows very rapidly from 0.1 M 0 to 3 during only 50,000 
yr. After Embryo 1 reaches the outer edge, it turns its migration direction 
inward. However, it encounters with another smaller embryo (Embryo 2) 
and inward migration is halted. There, Embryo 1 opens up a gap and starts 
slow outward migration due to distant encounters with inner planetesimals. 
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Embryo 1 is promoted to a full-embryo at t = 0.2 Myr during its inward 
migration. 

Embryo 2 and another embryo (Embryo 3) also migrate outward during 
the ontward migration of Embryo 1, and they also grow more rapidly than 
inner embryos, but not as rapidly as Embryo 1. They eventnally change their 
migration directions inward and come back to near the inner edge. Embryo 

2 occasionally attempts to migrate outward even after that, although its 
outward migration is not sustained. Two ontward migrations at t = 0.55 and 
0.65 Myr are halted by enconnters with another embryo at ~ 7 AU, which 
merges with Embryo 2 at t = 0.65 Myr. Embryo 2 is promoted to a fnll- 
embryo at this time. We do not clearly understand what prevents Embryo 2 
from snbseqnent outward migrations. The radial prohle of the surface density 
of planetesimals is no longer smooth due to occasional gap formation by 
embryos, and the outward migrations of Embryo 2 seem to be halted at local 
low density regions. Another possibility is that migrations of fnll-embryos are 
halted by stochastic effects since they are not massive enough as compared 
with tracers. This possibility may be tested using different transition masses 
from a sub-embryo to a full embryo. The masses of Embryo 2 and Embryo 

3 respectively reach 3.6 and 2.5 at 2 Myr. Their growth is slow at 
this time as they open np a gap in the planetesimal disk. 

Obviously, we need to perform more simulations with different parameter 
sets and conduct more careful analysis. Nevertheless, the example we showed 
here is a good demonstration that onr code can prodnce rapid growth of 
embryos during their outward migrations. 

Ngo (2012) □ performed simulations of accretion of cores of jovian planets 
using the LIPAD code. The initial conditions of his simulations are similar to 
ours. In his simulations withont fragmentation (his Figs. 5.14 and 5.15), sub¬ 
embryos are almost radially stationary even though the routine of sub-embryo 
migration is inclnded. Once sub-embryos are promoted to full-embryos, they 
suddenly start rapid radial migrations. It is nnclear to us why such large 
differences between snb-embryos and fnll-embryos are seen in his simulations. 
In spite of the large difference of embryo migration between his and onr 
simulations, the hnal masses of the largest embryos are close to each other. 


^His master thesis (Ngo, 2012) together with mpeg animations is available at 
https: / /qspace. library, queensu. ca/handle /1974/7400. 
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4. Discussion 

4 . 1 . Computational load 

The cpu time T^pu per timestep of orbital integration as a function of the 
number of embryos and the number of tracers Nt,. is roughly given as 



(65) 


where Ckep, (Tgrav, and Csta are the cpu times for a single orbital integration 
(a Kepler drift), for a single direct gravity calculation of a pair in the A^-body 
routine, and for a single set of calculations for a pair in the statistical routine, 
and iVint is the typical number of interlopers in the neighboring search region 
(Fig. 2). For a standard latest computer, Ckep ~ 4 x 10“^ s. From some of 
simulations, we hnd that C'grav/C'kep 0.15 Rlici C^sta/^kep ("Nj 0.3. 

As seen in Eq. ([H5D, the load on the statistical routine is controllable by 
modifying Ajnt and At/dt. Since the radial width of the neighboring search 
region is given by Eq. (jl]), oc for a fixed surface density. If we 

choose an azimuthal width as 66 oc becomes independent of Atr 

and the cpu time for the statistical routine can be linearly scaled with Ntr- 
Practically, one may choose a 66 so that the load on the statistical routine 
does not exceed the load on the A-body routine. 

The actual computational load on the statistical routine in some of our 
simulations is as follows. Recall that in the present paper, we call the statis¬ 
tical routine every 30 steps of orbital integrations {At = 306t). If the time 
step of the statistical routine is halved, its load is doubled. In the hybrid 
simulations shown in Fig. 5 (Atr = 200 and Aint ~ 10), the load on the statis¬ 
tical routine is only 10 % relative to the total load. The rest of 90 % is on the 
A-body routine and more than half of it is occupied by the Kepler equation 
solver. There is no direct gravity calculation in these simulations, as no em¬ 
bryo is included. In the hybrid simulation shown in Fig. 15 (Atr = 10000 and 
Aint ~ 100), the load on the statistical routine is about 60 % before embryos 
appear. If the number of embryos is larger than about 10, the load on the 
direct gravity calculation exceeds that for the Kepler equation solver, and 
the relative load on the statistical routine decreases. At 1.6 Myr of Fig. 15, 
at which 21 embryos exist, the load on the statistical routine is about 45 
%, that on the direct gravity calculation is 35 %, and that on the orbital 
integration is 20 %. Overall, the load of the statistical routine is comparable 
to or less than that on the A-body routine. 
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We also comment on parallelization. The current version of our code is 
serial. However, the direct iV-body routine is already parallelized and we test 
how much simulations are accelerated by parallelization in case without the 
statistical routine. In the tests, we £x the number of embryos to be 10 and 
vary the number of tracers. We find that simulations are accelerated only 
if the number of tracers is more than ~ 10^. In contrast, the pure iV-body 
simulations using the tree code, pkdgrav2, are accelerated by parallelization 
if N > 10^ (Morishima et al., 2010). The difference comes from the cpu 
times per step. To be benehtted by parallelization, the cpu time needs to be 
somewhat long to ignore latency. The computational cost per step for hybrid 
simulations is as small as those for iV-body simulations without particle- 
particle interactions. That ironically makes parallelization inefficient. 

4-2. Limitations of the code and future work 

We consider only random walk of tracers in tracer-tracer interactions. 
Thus, self-sustained migration occurs only to embryos in our code although 
small bodies (< rrito) may migrate in reality. To examine if migrations of 
small bodies play an important role for overall accretion history of embryos, 
high resolution simulations (with smaller rrito) are necessary. Minton and 
Levison (2014) showed that self-sustained migration occurs only if the em¬ 
bryo is much more massive than surrounding planetesimals and sufficiently 
separated from other embryos. Such conditions are usually fulhlled near the 
transition from the runaway growth stage to the oligarchic growth stage. The 
mass of the largest embryo at the transition is roughly given by 
(Ormel et al., 2010b; Morishima et ah, 2013), where ruiso is the isolation 
mass and mo is the initial mass of planetesimals. Therefore, if mo is very 
small (e.g., Weidenschilling, 2011), a very large number of tracers is nec¬ 
essary to represent an embryo at the transition by a single particle in our 
code. However, Minton and Levison (2014) also pointed out that the growth 
time scale needs to be longer than the migration time for an embryo to avoid 
encounters with other embryos during migration. If mo is very small, this 
criterion is not generally fulhlled for the largest embryo at the transition due 
to a very short growth time scale. Thus, masses of candidate embryos for 
self-sustained migration are probably not too small to be practically resolved 
by our code even if mo is very small. The argument here, of course, depends 
on disk parameters and radial location and needs more careful estimates. 

Our code can probably handle eccentric/inclined ringlets that are usually 
produced by external massive bodies (see Section 2.3). However, we have not 
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tested such simulations in this paper, and they remain for future work. A 
potential issue of our code in handling eccentric/inclined ringlets is accuracy 
of the global gravity forces (Section 2.7). Only for this calculation, we assume 
that a ring is axisymmetric and symmetric with respect to the invariant plane. 
These assumptions may not be appropriate if the radial/vertical thickness or 
the ringlet is much smaller than the radial/vertical excursion of individual 
particles. This is the case of some of narrow ringlets seen around solar system 
giant planets. If the global self-gravity plays an important role in structure 
evolution of narrow ringlets, we need a more advanced method to calculate 
it. 

The current code includes an option of hit-and-run bouncing (Appendix E) 
but not fragmentation. We plan to implement collisional fragmentation to 
our code in subsequent work. Since gas drag damps e and i of small frag¬ 
ments effectively, growth of embryos is accelerated by fragmentation of plan- 
etesimals, particularly in the oligarchic growth stage. On the other hand, 
since the total solid mass available for planetary accretion decreases due to 
rapid inward migration of fragments in the gaseous disk, the hnal masses 
of planets also decrease, unless planetary atmospheres effectively capture 
fragments. These conclusions were derived by simulations using statistical 
accretion codes (e.g., Inaba et ah, 2003; Chambers, 2008; Kobayashi et ah, 
2010, 2011; Amaro-Seoane et ah, 2014). In these studies, the torques on 
embryos exerted by small fragments were not taken into account. Levison et 
al. (2010) showed that small fragments shepherded by embryos push bash 
the embryos toward the central star, leading to rapid inward migration of 
the embryos. On the other hand, Kobayashi et al. (2010) pointed out that 
shepherded fragments are quickly ground down due to mutual collisions and 
very small ground fragments collide with or migrate past the embryos before 
pushing back the embryos. However, this has not been conhrmed by direct 
Lagrangian type simulations and is of interest in future work. 

5. Summary 

We developed a new particle-based hybrid code for planetary accretion. 
The code uses an A^-body routine for interactions with embryos while it can 
handle a large number of planetesimals using a super-particle approximation, 
in which a large number of small planetesimals are represented by a small 
number of tracers. Tracer-tracer interactions are handled by a statistical 
routine. If the embryo mass is similar to tracer masses and if embryo-tracer 


59 



interactions are handled by the direct iV-body routine, the embryo suffers 
artihcially strong kicks from tracers. To avoid this issue, sub-embryos are 
introduced. Accelerations of sub-embryos due to gravitational interactions 
with tracers are handled by the statistical routine whereas accelerations of 
tracers due to gravitational interactions with sub-embryos are handled by 
the direct A^-body routine. 

Our statistical routine hrst calculates the surface number densities and 
the orbital elements of interloping planetesimals around each target tracer 
(Section 2.2). Using the phase-averaged collision probability, whether a colli¬ 
sion between the interloper and the target occurs is determined (Section 2.4). 
If a collision occurs, the velocity changes are accurately calculated by match¬ 
ing the positions of the interloper and the target without changing their 
eccentricity and inclination vectors. Using the phase-averaged stirring and 
dynamical friction rates, the change rates of the orbital elements of the tracers 
are calculated (Section 2.5). These rates are converted into the accelerations 
of the tracers using both hrst and second order Gauss planetary equations. 
Planetesimal-driven migration of sub-embryos is basically handled by the 
A^-body routine but the migration rate is limited to the theoretical predic¬ 
tion derived by the statistical routine (Section 2.6). This unique routine can 
reproduce smooth, long-distance migrations of sub-embryos. 

We performed various tests using the new hybrid code: velocity evolution 
due to collisions, and that due to gravitational stirring and dynamical fric¬ 
tion, self-sustained migration of sub-embryos, formation of terrestrial planets 
both in the presence and absence of the gaseous disk, and formation of cores 
of jovian planets. All the test simulations showed good agreements with an¬ 
alytic predictions and/or pure A^-body simulations for all cases, except that 
the last test did not have a robust benchmark. The computational load on 
the portion of the statistical routine is comparable to or less than that for 
the A^-body routine. The current code includes an option of hit-and-run 
bouncing but not fragmentation that remains for future work. 
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Appendix A: Correction due to self-encounters 

In Section 2.2, we only consider encounters with interlopers. However, 
planetesimals in the target i may encounter with other planetesimals in the 
same target if ki > 1. To take into account this effect, we modify apparent 
planetesimal numbers in interlopers, if masses of interloping planetesimals 
are similar to those in the target tracer. For this purpose only, we assume 
axisymmetric distribution of planetesimals and introduce the mass bin for 
planetesimals. We also use the radial bin used in the neighboring search. We 
hrst calculate the total mass of tracers, radial and mass 

bins where the target locates. The summation is done for tracers in any 
azimuthal locations, excluding the target. Since the actual total mass seen 
from an individual planetesimal in the target tracer is + iki — ^)Tni, 

we use a corrected kj for the interloper. 


K = 


1 + 


{kj - l)mi 
Ef'kjTTlj 


Vj, 


( 66 ) 


if it is in the same radial and mass bins with the target. If there is no other 
tracer in the same radial and mass bins, the correction is not applied. We 
hnd that this self-encounter effect is very small at all as we do not identify 
any difference in simulations with and without the correction. 

The effect can also be checked by simulations with a same total number of 
planetesimals but using different numbers of tracers. Although a simulation 
with a small number of tracers has larger statistical fluctuations in various 
physical quantities, such as the mean eccentricity, any differences in time- 
averaged quantities are not well identihed. Thus, we can omit this correction, 
and mass bins are unnecessary in our code. 


Appendix B: Non-dimensional collision, stirring, and migration 
rates 

With Hill’s approximation, the equations of motion for the relative mo¬ 
tion between the target and the interloper are reduced to Hill’s equations 
(Hill 1878; Nakazawa et ah, 1989). The time and the distance of Hill’s equa¬ 
tions can be normalized by the inverse of orbital frequency, and aijhij, 
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respectively, where aij = {ai + aj)/2 is the mean semimajor axis and hij is 
the reduced mutual Hill radius [Eq. ([H])]. The solution of the relative motion 
in the absence of mutual gravity is expressed by non-dimensional relative 
orbital elements, b, e^, ir, rur and Hr- The elements b, and h are given as 


b 


^ij hij 





(67) 


where dij = aj — ai is the difference of the semimajor axes, the relative 
eccentricity Cr and the relative inclination v are dehned in Eq. (E]) together 
with the phases cUr and Hr. 


B.l. Collision rates 

The non-dimensional collision rate Pcoi{c,ir) is defined as 

/~~\ /* ^111 doj rdhlr / \ 

T’coi(er,v) = J Pco\-\b\db , ( 68 ) 

where Pcoi is unity for collisional orbits and zero otherwise. 

In the dispersion dominant regime [(e^ -|- > 2], Pcoi is given as 

(Greenzweig and Lissauer, 1990; Inaba et ah, 2001) 


Tcol,high(Cr, ^r) — 


2rr 


e? it 


TT 


^(0 + 


6 K{C) 
bp e? + P 


(69) 


where hp = {si + Sj)/{aijhij) with the physical radii Si and Sj, and E{Q 
are the complete elliptic integrals of the hrst and second kinds, respectively, 
and CY = 3e^/(4e^ 4i^). 

In the shear dominant regime [(e^ + < 2], if the inclination is small 

enough (i-^ < bp^^^; see Goldreich et ah (2004) for this criterion) all pairs of 
particles entering the mutual Hill sphere collide, and Pcoi becomes indepen¬ 
dent of hr and h as (Ida and Nakazawa, 1989) 


Tcol,low(€'rGr) 11.3 -y^T^p. (”^0) 

In the shear dominant regime with moderate inclinations (bp^^^ < i < 2), 
Pcoi depends only on i as (Ida and Nakazawa, 1989) 

~ 14Tt* 

Tcol,med(brGr) ~ 

47r«r 
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Following Inaba et al. (2001), we connect the collision rates in different 
regimes as 


^col(er, ir) = min[Pcol,med, (^"col^high + ^col^low) (^2) 

The collision rate given by this expression well agrees with that derived by 
numerical integrations (Ida and Nakazawa, 1989; Greenzweig and Lissauer, 
1990). 


B.2. Stirring rates 

The non-dimensional stirring and dynamical friction rates are defined as 


h) 

QyS (^d h) 
-Pdf(G, h) 



(73) 

(74) 

(75) 


where 57^, Si^, and Ssj. are the changes of e^, and Bj. [Eq. ([6])] during a 
single encounter. The non-dimensional dynamical friction rate for inclination 
is obtained by replacing (dr • Se^/e^) by (ir • Sir/i^) in Eq. fl75ll and is known 
to be the same with Pdf (Tanaka and Ida, 1996). 

The rates for the shear dominant regime are (Ida 1990, Ohtsuki et ah, 

2002 ) 


A^S,low — 

75.4, 

(76) 

QvS,low = 

2il -I- 5efir, 

(77) 

TdFjIow ~ 

10. 

(78) 


For Qvs.iow, we adopt a functional form that is similar to the one with the 
Rayleigh distribution (Ohtsuki et ah, 2002) but the factors are modified. 
The rates for the dispersion dominant regime are given as (Tanaka and 
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Ida, 1996) 


-PvS.high — 


QvS.high — 


36 


P^ 


DFjhigh 


5^(0 - 3r^^(C) 


■nir\Je^ + ir ^ 
X In (A^ + 1), 
36 


+ 4*2 


12*2 

K{Q - 


'Kir\J + *2 L 

X In (A^ + 1), 


e2 + 4*2 


7r*i.(e2 + Ai‘^)\ + *2 


-.E{0 ln(A2 + l), 


(79) 


(80) 

(81) 


where 

A = v(e1 + *?)/3. (82) 

In Pvs.high and Qvs,high, we omit the term Is?/{Is? + 1) as it is negligible 
compared with In (A2 + 1) in the dispersion dominant regime. 

We connect the stirring and dynamical friction rate at the shear and 
dispersion dominant regimes in a similar manner of Ohtsuki et ah (2002): 


The coefficients are 


Pvs 

— Oi Pvs, low + Pvs,high 

(83) 

Qys 

= 02(5vS,1ow + QvS,high 

(84) 

Pdf 

~ C*3pDF,low T Pdf, high- 

(85) 

Cl = 

ln(A2 + 1)/A2, 

(86) 

c, = 

ln(10A2er + l)/(10A2er), 

(87) 

Cs = 

ln(10A2 + l)/(10A2). 

(88) 


The approximated formulae well agree with the rates derived by three-body 
orbital integrations [Fig. 8 of Ida (1990) and Fig. 1 of Ohtsuki et ah (2002)], 
except for the region with Cr > 4 and *r < 1 (see Section 3.2). 


B.3. Migration rates 

The non-dimensional migration rate Rmg for one-sided torque is dehned 


as 



dev 

-do, 


(27r) 


(89) 
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where Sb is the change of b during a single encounter and the negative sign 
represents the back reaction on the embryo. At = 0, Rug = 9.4 [Ida 

et ah, 2000; their Eq. (10)]. It is expected that Rug in the shear dominant 
regime is about the same value as the migration rate is almost constant in 
this regime (Kirsh et ah, 2009): 


.Rmg,1ow — 9.4. 


(90) 


For the dispersion dominant regime, Ida et al. (2000) derived the migration 
rate averaged over and fir as [their Eq. (18)] 


dzuT-dfli 

(27r)2 


6 b 


12 b In(A2 + l) 

1^1 ^(e?-&2)(e2+ •2_3^2/4)3 


Inserting this equation into Eq. fl8^ . we have 


(91) 


R 


MG,high 


72 In (A^ + 1) gg 
’"tSr /g2 + |2 e'? + 4!?’ 


(92) 


where the integration range of b is limited between 0 and Cr, and In (A^ + 1) 
is approximately treated as a constant in the integration. Analogous to the 
stirring rates, we connect the migration rates at the shear and dispersion 
dominant regimes but the value is assumed not to exceed i?MG,iow 


RmG — MIN[(Cii?MG,low + .RMG,high), .RmG,1ow]- (93) 

While this expression must be correct in the low and high velocity limits, 
its accuracy at the intermediate regime needs to be confirmed by comparing 
with the rate derived from three-body orbital integrations. 

Appendix C: Splitting of dynamical friction 

When the non-dimensional collisional and stirring rates are derived, uni¬ 
formities of zUr and fir are assumed. However, we do not assume uniformities 
of the relative longitudes, tUij = wj — Wi and 12^ = Vtj — Hj. Alignments 
of the longitudes occur, for example, if the system is perturbed by a giant 
planet and the forced eccentricity is much larger than the free eccentricities. 
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Ida (1990) showed that the eccentricity evolution of a particle i due to 
dynamical friction of particles j is 



where the characters appeared in these equations are dehned around Eqs. (flSD 
and If Wij's of interlopers are randomly distributed and the distribu¬ 

tion of Cj is the Rayleigh distribution, the term with cos Wij becomes small 
enough after averaging over interlopers (Ida, 1990). The equilibrium state for 
this case leads to the well-known energy equipartitioning with 

dynamical friction only. With viscous stirring, however, the energy equipar¬ 
titioning is achieved only if the mass distribution is very steep (Rahkov, 
2003). 

If the distribution of Wij is not random, as may be the case of a forced 
system, the term with cos Wij in Eq. (1951) cannot be ignored. In a forced 
system, the eccentricity vectors may be split as e* = ef+Bio and Bj = Bf + Bjo, 
where Cf is the forced eccentricity vector, Cjo and b^q are the free eccentricity 
vectors, respectively (Murray and Dermott, 1999). Inserting these forms 
into Eq. fl9T)) and setting ef • (e^o — e^o) = 0 due to the uniformity of ciJr, we 
can End an equation similar to Eq. fl9^ but e* and Cj replaced by the free 
eccentricities ejo(= |ejo|) and ejo(= |ejo|) and Wij replaced by the relative 
longitude of the free eccentricities. Since there is no correlation between the 
directions of e^o and Cjo, the third term in the square bracket of the r.h.s of 
Eq. fl9S]) can be neglected. Thus, dynamical friction tends to lead the energy 
equipartitioning of the free eccentricities, if Cf is shared for all particles. A 
similar conclusion can be derived for inclinations. 

Now we consider splitting Eq. fl9^ into two terms: a stirring term that is 
preferably larger than the viscous stirring term [as required in Eq. 0341) ] and 
a damping term that decreases the free eccentricity. We hnd the forms that 
can fulhll such conditions as 



(96) 


(97) 
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The same equations are given in Eqs. (1T5]) and ffT^ . In a forced system, 
the terms in the parenthesis of the r.h.s. of Eq. (in7|) are reduced to e^Q 
using Bf ■ (sjo — Bio) = 0 and Cjo ■ e^o = 0 due to the uniformity of cUr. 
The l.h.s. of Eq. fl^ is also reduced to def^ldt because the orientation of 
Cjo changes randomly due to gravitational interactions with other tracers 
[Bf ■ (Bio/dt) = 0]. Thus, the damping term [Eq. ((971)] decreases the free 
eccentricity, as required. Usually in astrophysics, only this damping term is 
called dynamical friction (Binney and Tremaine, 1987). A similar splitting is 
possible for inclinations, and the split terms are given by Eqs. fl2T]) and ([22]). 


Appendix D: Radial diffusion coefficient 


The classical theory of random walk shows that the diffusion coefficient Di 
for the one dimensional problem is {(Aaj)^)/(2At), where Aa, is the change 
in the semimajor axis for the tracer i during the time step At and the angle 
brackets mean averaging over multiple changes. The change of a* due to a 
single encounter with other planetesimal is defined to be Sat = Ujttijhijdb, 
where 6 b is again the change of b during an encounter. Thus, the diffusion 
coefficient {Di)j for the tracer i due to encounters with planetesimals in the 
tracer j is given as 

(A), = (98) 

where Ab is the change of b during At and the subscript for the angle brackets 
means averaging over planetesimals in the tracer j. 

Petit and Henon (1987) performed three-body orbital integrations and 
showed that {b6b)/{{6bY) = —3.07/17.72 ~ —1/6 for Cr = h = 0. Although 
it is unclear if this factor is applicable for er,ir > 0, we use the relationship 
for any cases. This gives {{6bY) = 3{6b'^)/2, using {db^ = —2b6b + 5b'^. The 
averaged change ((A6)^)j is given by 


((A&)^), 

At 


2 u 2 

Uja^jhij uk 




dj caj 

7 ^ 


212 ^0^^ hi j‘^UK{PvS + Qvs), 


(99) 


where we used 56^ = (4/3)(5e^ + 6i‘^) that is derived from the Jacobi integral 
of Hill’s equations [Eq. fH3|) ]. and Pys and Qvs are dehned in Eqs. fl73l) and 
dZH). 
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Inserting Eq. (1^ into Eq. the diffusion coefficient {Di)j is given as 

{Di)j = i/j afjhfjUKiPvs + Qvs)- ( 100 ) 

Finally, the diffusion coefficient due to encounters with all interlopers is given 
by Di = Ohtsuki and Tanaka (2003) derived the viscosity for the 

equal-mass planetesimals disk and their viscosity is (2/9)D*. The diffusion 
coefficients derived by Ormel et ah ( 2012 ) and Glaschke et ah (2014) are 
similar to ours. 

Appendix E: Hit-and-run collisions 

In Section 2.4, we show how the total number of tracers changes due to 
merging between planetesimals. On the other hand, in hit-and-run collisions, 
there should be no change in the number. Before describing how to handle 
hit-and-run collisions between planetesimals in two tracers, we hrst show how 
to handle a hit-and-run collision between embryos or between a full embryo 
and a tracer. 

E.l. Criterion of hit-and-run collisions 

Consider an impact between the target with mass rui and the impactor 
with mass mj (< mi). Let the impact velocity and the impact angle be Uimp 
and 6 c {6c = 0 for a head-on collision). If Uimp > Ucr, where Vcr is the critical 
velocity, the collision results in escape of the projectile after bouncing on the 
target. Based on thousands of SPH simulations, Genda et ah (2012) derived 
the formula of Vcr as 

^ = cir©'^^ + caT + C 30 '=® + C4, (101) 

Vesc 

where Uesc is the mutual escape velocity, T = {mi — mj)/{mi -\- mj), and 
0 = 1— sin 6 *c. The coefficients are Ci = 2.43, C2 = —0.0408, C3 = 1.86, 
C 4 = 1.08, and C 5 = 5/2. 

If the impact is judged as a hit-and-run collision (uimp > Vcr), the post¬ 
impact velocities of the projectile and the target are determined as follows. 
We do not consider any mass exchanges between the target and the impactor 
for simplicity following Kokubo and Genda (2010). The impact velocity vec¬ 
tor is decomposed to the normal and tangential components {vn and Vt) 
relative to the vector pointing toward the center of the target from the cen¬ 
ter of the impactor. The tangential component of the post-impact relative 


velocity v[ is assumed to be the same as Vt- The normal component of the 
post-impact relative velocity is given as 



(for v[ > Uesc), 

(otherwise). 


( 102 ) 


The velocity change described here is similar to but slightly different from 
that of Kokubo and Genda (2010), who always set v'^ = 0 and adjust 

The normal and restitution coefficients, €„ and et, are defined as €„ = 
—v'^/v-a and Ct = v[/vt- While the above case basically assumes Cn = 0 and 
Ct = 1, these coefficients can be changed to any values in simulations (in 
Section 3.1, we use €„ = 0 and 0.8). 

E.2. Procedure of hit-and-run collisions between tracers 

Now we consider hit-and-run collisions between planetesimals in two trac¬ 
ers. We will handle them as consistent as possible with those between em¬ 
bryos. If a collision is judged to occur in the statistical routine, the target and 
the impactor are moved to the impact position without changing their ec¬ 
centricity and inclination vectors, as described in Section 2.4.2. The relative 
velocity, = Ini.!, at inhnity (without acceleration due to mutual gravity) 
is calculated at this impact position. If kjb > 0, only some fraction of plan¬ 
etesimals in the tracer j are involved in the collision (Fig. 4). The impact 
velocity is given as Uimp = (n^ -|- We ignore deflection of the relative 

orbit during acceleration due to mutual gravity so that the direction of the 
impact velocity vector is identical to that of the relative velocity vector v^. 
The problem is that we cannot directly use the relative position vector of 
the tracers because their mutual separation is not exactly the sum of the 
radii of a colliding pair after matching their positions. Instead of using the 
actual relative position vector, we make a fictitious relative position vector 
rrf. The hctitious impact angle 6 c is given by cos^c = (ur ■ rrf)/(nr|ri.f|). 
Assuming there is no correlation between and r^f, the probability that the 
impact angle is between 6 c and 6 c + d6c is given as df = 2 sin 6 c cos 6 cd6c (e.g.. 
Shoemaker and Wolfe, 1982). We choose an arbitrary 6 c between 0 and 90 
degrees following this distribution and also choose the tangential direction of 
Vj-f relative to at random. Then, we judge whether the collision results in 
merging or bouncing using Eq. fllOll) in which mj is replaced by Am*. If it 
is merging, we follow the process described in Section 2.4.2, using v,. as the 
impact velocity vector. 
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If it turns out to be a hit-and-run collision, we calculate the post impact 
velocities using v[ = Vt and v'^ from Eq. fll02p . Then their velocities are 
modihed so that their relative velocity at inhnity is given as v[ = {v[^ + v'^ — 
again ignoring deflection of the orbits due to mutual gravity. If the 
interloper is split into two tracers [Case (a) in Eq. flT^ ] before the collision, 
we merge them by summing their momenta after matching their positions. 
The merged tracer has the same rrtj and kj as those of the original interloper 
before splitting. The semimajor axis of the interloper is adjusted so that 
the z-component of the total angular momentum of is conserved. Einally, 
the mean longitudes of the target and interloper are placed to the original 
values. 

Appendix F. Duration of a tracer in the neighboring search region 
around a sub-embryo 

Here we explain how to derive the correction factor Tj/T^ in appeared in 
Eq. dlH]). Let ns dehne the inner and outer radii of the region i for neighboring 
search around the embryo be rjn and rout (Fig. 2). The interloper j has the 
semimajor axis a^, the eccentricity e^, and the inclination ij. The maximum 
and minimum distances of the interloper j, rmax and rmin, from the central 
star projected on the invariant plane {z = 0) are given by 

’’max = a j{l + ej) cos ij (103) 

^max = aj{l- Cj) cosij. (104) 

If ’’max > ’’out, we dehne the eccentric anomaly Eont (< tt) at r = rout as 

’’out = aj{l- Cj cos Eont) cosij (105) 

and the mean anomaly is also dehned as Mout = .^out — Cj sinEout- If ’’max < 
’’out, we set Mout = ”. In a similar manner, if rmin < ’’in, we dehne the 
eccentric anomaly E^n (< ’’) at r = r^ as 

’’in = 0^(1 — Cj cos Fin) COS Zj (106) 

and the mean anomaly Mjn = Ein — Cj sinEju. If rmin > ’’in, we set M^ = 0. 
Then, the period Tj in in which the interloper j is radially in the region i 
during the orbital period Tj is given as 






Appendix G: Comparison with Levison et al. (2012) 

Here we summarize the differences between the LIPAD code (LDT12) 
and our new code. 

Spatial grid for neighboring search 

In the LIPAD code, to derive the spatial density of planetesimals, tracers 
are placed in mass and radial grids and the vertical distribution is assumed to 
be Gaussian. Our code has neither these hxed grids nor the assumption for 
the vertical distribution (except for the global gravity calculation described 
in Section 2.7). Thus, our code has Lagrangian natures more than the LIPAD 
code, although these differences are probably not fundamental for the test 
simulations shown in this paper. Our scheme is more benehcial for systems 
with non-axisymmetric and/or inclined structures. 

The advantage of the assumption of the vertical distribution in the LIPAD 
code is that changes in collision and stirring rates at different are explicitly 
taken into account. Changes in collision rates with z are taken into account 
in our model too, although implicitly, as we move both the target and the 
interloper to the collisional point (Section 2.4); collisions occur more often 
near the mid-plane than at high \z\. In the stirring routine (Section 2.5), our 
method ignores this effect and this probably introduces some inaccuracies in 
orbital evolutions of individual particles. Nevertheless, the velocity evolution 
of the system as a whole is reproduced reasonably well by our method as 
shown in the test simulations (see more discussion for viscous stirring and 
dynamical friction given below). 

Number of planetesimals in a tracer 

The number of planetesimals in a tracer, /cj, is an integer in our model 
(Section 2.1), contrary to LDT12. If a tracer with a non-integer number of 
planetesimals is promoted to a single sub-embryo, the extra mass, {ki — l)mj 
(if 1 < ki < 2), needs to be incorporated into the sub-embryo or transferred 
to other tracers which may have very different planetesimal masses. Both are 
not accurate although it is not clear to us how it is handled in the LIPAD 
code. 

Orbital isolation 

In our code, orbital isolation between planetesimals does not occur as we 
choose a narrow radial half width Sr of the neighboring search region given 
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by Eq. (|1]). If we adopt a larger 6 r, orbital isolation can potentially occur. 
In this case we may need to consider how isolated bodies interact through 
distant perturbations. In general, however, the radial resolution is sacrihced 
with a large 6 r. In our current scheme, interactions between runaway bodies 
are approximately represented by occasional close encounters rather than 
distant perturbations. 

The LIPAD code explicitly handles orbital isolation and distant pertur¬ 
bations (e.g., their Fig. 10). This means that they use the radial grid size 
larger than ours. Even with a large Sr, the exitisting Lagrangian codes can¬ 
not handle orbital isolation of planetesimals that are much less massive than 
rritQ, because a small number of runaway bodies cannot be represented by 
tracers with nearly hxed masses. 

Collisions between planetesimals and sub-embryos 

Collisions between tracers and sub-embryos are handled in the A^-body 
routine in the LIPAD code. Thus, after the promotion the hrst impact of 
a tracer doubles the mass of a sub-embryo. This makes artihcial kinks in 
mass spectra at the transition mass, m^o- The influence of this effect on 
evolution of the overall system is small as discussed by LDT12, particularly 
if the masses of the largest embryos are much larger than mto- 

In our method, collisions of planetesimals with sub-embryos are handled 
by the statistical routine (Secs. 2.4 and 2.6). The mass spectra are smooth at 
mtQ (see Fig. 12), because only some of planetesimals in a tracer collide with 
a sub-embryo. Note that we assume uniformities in phase space (Section 2.3) 
that seem to be valid at least in test simulations shown in this paper, while 
no assumption for phase-space is necessary for the A^-body routine. 

Collision velocity between tracers 

The impact velocity between tracers is calculated more accurately in our 
code than the LIPAD code. As described in Section 2.4, we estimate the 
impact velocity between tracers by matching their positions without changing 
their eccentricity and inclination vectors. 

In the LIPAD code, only the interloper’s longitude is aligned to that of 
the target by rotating the interloper’s coordinates, not letting the interloper 
move along its Keplerian orbit. This changes the longitude of pericenter 
of the interloper by the rotation angle. LDT12 managed to minimize this 
issue by choosing the azimuthally closest interloper for the collision with the 
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target, as originally done in Levison and Morbidelli (2007). On the other 
hand, our method allows any interlopers to collide with the target. 

The timescale of collisional damping shown in Fig. 2 of LDT12 is clearly 
too short, probably because they described the initial condition of a wrong 
simulation. 

Viscous stirring and dynamical friction 

The largest difference between our code and the LIPAD code probably 
appears in the methods to handle viscous stirring and dynamical friction. In 
the LIPAD code, the probability that an interloper enters within a mutual 
Hill radius of the target is calculated. If an encounter occurs, the acceleration 
of a smaller body of the pair is given by solving a pair-wise gravitational 
interaction. The acceleration of a lager body of the pair is not calculated 
in this routine, but the damping forces are given by the Chandrasekhar’s 
formula of dynamical friction. Instead of taking into account stirring on the 
larger body, they set lower limits of e and i of the larger body based on 
energy equipartitioning. Some inaccuracies seem to be introduced in their 
approach. First, energy equipartitioning is not usually achieved except with 
a very steep mass spectrum (Rahkov, 2003). Second, the Chandrasekhar’s 
formula of dynamical friction is not applicable to the shear-dominant regime 
(e.g., Ida, 1990; Ohtsuki et ah, 2002). 

In our method, we do not solve individual pair-wise gravitational in¬ 
teractions. Instead, we calculate the time averaged changes of the orbital 
elements due to multiple encounters with tracers, using the phase-averaged 
stirring rates (Section 2.5.1). These changes of orbital elements are then con¬ 
verted into the accelerations of tracers (Section 2.5.2). This approach allows 
us to handle viscous stirring, the stirring and damping parts of dynamical 
friction due to interactions between all pairs in the same manners regardless 
of their mass ratios. Our approach is based on uniformities of the phases 
(Section 2.3), and methods which explicitly solves pair-wise gravitational in¬ 
teractions (i.e., the acceleration of a smaller body in the LIPAD code) are 
probably more accurate than ours. We would like to point out, however, 
that the LIPAD code also partly assumes uniformities of the phases because 
it uses the phase-averaged collision/encounter probability of Greenzweig and 
Lissauer (1990); note that Pcoi we use and Fg in the LIPAD code are the 
same quantity in different formats. 
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Planetesimal-driven migration of sub-embryos 

LDT12 included planetesimal-driven migration of sub-embryos by adding 
torques on them calculated by a series of three-body integrations (Sun, the 
target, and the interloper) during LIPAD simulations. Each interloper is 
randomly chosen from one in seven Hill radii of the target. It has the same 
semimajor axis, orbital eccentricity, and inclination as one of the tracers but 
its phases are modified. As far as we understand, their point is that the torque 
on the sub-embryo is smoothed by phase averaging. This may be true in 
the dispersion dominant regime {eij/hij,iij/hij 3> 1). In the shear dominant 
regime, however, the phase dependence of the torque is weak and the averaged 
toque turns out to be as strong as the torque from a single interloping tracer 
at all [see Fig. 8 of Ohtsuki and Tanaka (2003), for example]. Thus, sub¬ 
embryos should occasionally suffer strong kicks from tracers. Although it is 
unclear to us whether the point described above is indeed an issue for the 
LIPAD code, sub-embryos in LIPAD simulations (Ngo, 2012) look radially 
much more stationary than those in our simulations (Section 3.6). 

In our method, the angular momentum change of a sub-embryo is stored 
in the A^-body routine but its release rate is limited to the theoretical pre¬ 
diction derived by the statistical routine (Section 2.6.2). Theoretical justih- 
cations of our method are discussed in Section 3.3. 

Three-body integrations 

The LIPAD code performs a series of three-body integrations during a 
simulation for two cases. The first one is for the routine of sub-embryo mi¬ 
gration. The another one is for the stirring routine in the shear dominant 
regime. Our code does not need to perform three body integrations. Migra¬ 
tion of sub-embryos are handled as explained above. Stirring in the shear 
dominant regime is handled using the approximated formulae of the phase- 
averaged stirring rates (Section 2.5 and Appendix B). 
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